This report describes the results of a preregistered study available at: https://osf.io/gkd8s.
Note also that this data has been cleaned beforehand. Several
datasets over three measurement times were merged (joined) through an
inner join so as to keep only participants who at least participated at
each step of the study. Missing data will be imputed later on.
Duplicates were addressed with the rempsyc::best_duplicate
function, which keeps the duplicate with the least amount of missing
values, and in case of ties, takes the first occurrence. However, for
duplicate participation in the activities and exercises, rather than the
first occurrence, the occurrence with the higher completion percentage
(% of total activity time) was taken instead.
library(rempsyc)
library(dplyr)
library(interactions)
library(performance)
library(see)
library(report)
library(datawizard)
library(modelbased)
library(ggplot2)
library(bestNormalize)
library(psych)
library(visdat)
library(missForest)
library(doParallel)
library(ggplot2)
library(emmeans)
library(sjlabelled)
library(eefAnalytics)
library(tidyr)# Read data
data <- readRDS("Data/finaldataset_n381.rds")
report_participants(data, threshold = 1) %>% cat381 participants (Mean age = 24.9, SD = 4.1, range: [18, 35], 43.0% missing; Gender: 45.1% women, 11.8% men, 0.00% non-binary, 43.04% missing)
# Allocation ratio
report(data$T1_Group)x: 3 levels, namely Meditation (n = 65, 17.06%), Reflection (n = 56, 14.70%), Waitlist (n = 96, 25.20%) and missing (n = 164, 43.04%)
report(data$T2_Condition)x: 2 levels, namely Control (n = 119, 31.23%), Depleted (n = 98, 25.72%) and missing (n = 164, 43.04%)
In this section, we are preparing the data for analysis: (a) taking
care of preliminary exclusions, (b) checking for and exploring missing
values, (d) imputing missing data with missForest, (e)
computing scale means, and (f) extracting reliability indices for our
scales.
Here, we report on participant exclusions and corresponding changes on the sample.
First, we want to exclude those who did not get assigned a group. How many people would we exclude?
data %>%
filter(is.na(T1_Group)) %>%
nrow()## [1] 164
164 people. Investigation of these data suggests these participants did not complete the questionnaires and did not reach the group assignment. It is safe to exclude them, so let’s do so.
data <- data %>%
filter(!is.na(T1_Group))
report_participants(data, threshold = 1) %>% cat## 217 participants (Mean age = 24.9, SD = 4.1, range: [18, 35]; Gender: 79.3% women, 20.7% men, 0.00% non-binary)
# Allocation ratio
report(data$T1_Group)## x: 3 levels, namely Meditation (n = 65, 29.95%), Reflection (n = 56, 25.81%)
## and Waitlist (n = 96, 44.24%)
report(data$T2_Condition)## x: 2 levels, namely Control (n = 119, 54.84%) and Depleted (n = 98, 45.16%)
Second, we know that we only want to keep participants who had a participation level of at least 2/3 of all activities and exercises. Let’s see the distribution of participants’ participation, by group. However, we do not want to exclude participants from the control group, so we will give them an artificial 100% participation rate. How many people would we exclude?
data <- data %>%
mutate(part.percent = convert_na_to(part.percent, 1))
data %>%
filter(part.percent < 2/3) %>%
count(T1_Group)| T1_Group | n |
|---|---|
| Meditation | 6 |
| Reflection | 1 |
7 people. Let’s exclude them.
data2 <- data
data <- data %>%
filter(part.percent >= 2/3)
report_participants(data, threshold = 1) %>% cat## 210 participants (Mean age = 24.8, SD = 4.0, range: [18, 35]; Gender: 80.5% women, 19.5% men, 0.00% non-binary)
# Allocation ratio
report(data$T1_Group)## x: 3 levels, namely Meditation (n = 59, 28.10%), Reflection (n = 55, 26.19%)
## and Waitlist (n = 96, 45.71%)
report(data$T2_Condition)## x: 2 levels, namely Control (n = 113, 53.81%) and Depleted (n = 97, 46.19%)
Let’s also exclude those who failed 2 or more attention checks (i.e., keep with those with a score of two or more).
data <- data %>%
mutate(att_check = rowSums(
select(., T1_attention1, T2_attention1, T3_attention1), na.rm = TRUE))
data %>%
count(att_check)| att_check | n |
|---|---|
| 0 | 3 |
| 1 | 4 |
| 2 | 13 |
| 3 | 190 |
There’s 7 more exclusions here.
data <- data %>%
filter(att_check >= 2)Here, we report various sample demographics.
report_participants(data, threshold = 1) %>% cat203 participants (Mean age = 24.8, SD = 3.9, range: [18, 35]; Gender: 81.3% women, 18.7% men, 0.00% non-binary)
report_participants(data, threshold = 1, group = "T1_Group") %>% catFor the ‘T1_Group - Meditation’ group: 58 participants (Mean age = 25.2, SD = 4.5, range: [18, 35]), for the ‘T1_Group - Reflection’ group: 53 participants (Mean age = 25.1, SD = 3.7, range: [19, 35]) and for the ‘T1_Group - Waitlist’ group: 92 participants (Mean age = 24.4, SD = 3.7, range: [18, 35])
# Allocation ratio
report(data$T1_Group)x: 3 levels, namely Meditation (n = 58, 28.57%), Reflection (n = 53, 26.11%) and Waitlist (n = 92, 45.32%)
report(data$T2_Condition)x: 2 levels, namely Control (n = 113, 55.67%) and Depleted (n = 90, 44.33%)
“How did you hear about the study?”
get_label(data$T1_recruitment) %>% cat## Comment avez-vous entendu parler de l'étude? - Selected Choice
report(data$T1_recruitment)## x: 8 levels, namely Email (n = 49, 24.14%), Facebook (n = 116, 57.14%), google
## (n = 5, 2.46%), Kijiji (n = 8, 3.94%), Other (n = 0, 0.00%), SQRP bulletin
## board (n = 10, 4.93%), SSMU Marketplace (n = 6, 2.96%) and Word of mouth (n =
## 9, 4.43%)
data %>%
count(T1_recruitment, sort = TRUE)| T1_recruitment | n |
|---|---|
| 116 | |
| 49 | |
| SQRP bulletin board | 10 |
| Word of mouth | 9 |
| Kijiji | 8 |
| SSMU Marketplace | 6 |
| 5 |
“Have you already completed a psychology course?”
get_label(data$T1_psycho.class) %>% cat## Avez-vous déjà suivi un cours de psychologie?
data %>%
count(T1_psycho.class, sort = TRUE)| T1_psycho.class | n |
|---|---|
| Yes | 124 |
| No | 79 |
“Have you ever tried virtual reality?”
get_label(data$T1_virtual.reality) %>% cat## Avez-vous déjà essayé une sorte de réalité virtuelle?
data %>%
count(T1_virtual.reality, sort = TRUE)| T1_virtual.reality | n |
|---|---|
| No | 137 |
| Yes | 66 |
“What is your meditation experience?”
get_label(data$T1_medi.experience) %>% cat## Quelle expérience de méditation avez-vous?
data %>%
count(T1_medi.experience, sort = TRUE, .drop = FALSE)| T1_medi.experience | n |
|---|---|
| < 5 hours | 161 |
| Between 5 and 10 hours | 42 |
| > 10 hours | 0 |
“Do you have a history of psychiatric or neurological disorders?”
get_label(data$T1_disorders) %>% cat## Avez-vous des antécédents de troubles psychiatriques ou neurologiques?
data %>%
count(T1_disorders, sort = TRUE)| T1_disorders | n |
|---|---|
| No | 203 |
“Do you have normal or corrected-to-normal vision?”
get_label(data$T1_vision) %>% cat## Avez-vous une vision normale ou corrigée à la normale?
data %>%
count(T1_vision, sort = TRUE)| T1_vision | n |
|---|---|
| Yes | 203 |
“Do you own a smart phone of type iPhone or Android?”
get_label(data$T1_phone) %>% cat## Possédez-vous un téléphone intelligent de type iPhone ou Android?
data %>%
count(T1_phone, sort = TRUE)| T1_phone | n |
|---|---|
| Yes | 201 |
| No | 2 |
“Do you live in Quebec at this time?”
get_label(data$T1_quebec) %>% cat## Habitez-vous au Québec présentement?
data %>%
count(T1_quebec, sort = TRUE)| T1_quebec | n |
|---|---|
| NA | 168 |
| Yes | 34 |
| No | 1 |
get_label(data$T1_student) %>% cat## Êtes-vous étudiant.e?
data %>%
count(T1_student, sort = TRUE)| T1_student | n |
|---|---|
| Student | 160 |
| Non-student | 43 |
get_label(data$T1_student.program_cat) %>% cat## Dans quel programme étudiez-vous?
report(data$T1_student.program)## x: 111 entries, such as psychologie (11.33%); baccalauréat en psychologie
## (2.96%); criminologie (2.96%) and 108 others (43 missing)
data %>%
count(T1_student.program_cat, sort = TRUE) %>%
filter(n > 3)| T1_student.program_cat | n |
|---|---|
| Psychology | 45 |
| NA | 43 |
| Other | 29 |
| Biology | 7 |
| Teaching | 7 |
| Biomedical Sciences | 6 |
| Law | 5 |
| Nursing | 5 |
| Social Work | 5 |
| Arts | 4 |
| History | 4 |
| Neuroscience | 4 |
| Sexology | 4 |
get_label(data$T1_workplace) %>% cat## Dans quel domaine travaillez-vous?
report(data$T1_workplace)## x: 39 entries, such as informatique (1.48%); administration (0.99%);
## psychologie (0.99%) and 36 others (160 missing)
data %>%
count(T1_workplace_cat, sort = TRUE) %>%
filter(n > 1)| T1_workplace_cat | n |
|---|---|
| NA | 160 |
| Other | 12 |
| Health | 5 |
| Information and Technology | 5 |
| Psychology | 5 |
| Admin | 4 |
| Arts & Culture | 4 |
Participants were asked questions about their meditation practice during the 13 weeks.
get_label(data$T3_post.medipractice) %>% cat## Avez-vous pratiqué la méditation au cours des 13 dernières semaines?
data %>%
count(T3_post.medipractice, sort = TRUE)| T3_post.medipractice | n |
|---|---|
| Non | 107 |
| Oui | 95 |
| NA | 1 |
get_label(data$T3_medipractice.which) %>% cat## Quel type de méditation avez-vous pratiqué? - Selected Choice
report(data$T3_medipractice.which)## x: 21 entries, such as NA (53.20%); mindfulness (20.69%); LKM (7.88%) and 18
## others (0 missing)
data %>%
count(T3_medipractice.which, sort = TRUE)| T3_medipractice.which | n |
|---|---|
| NA | 108 |
| mindfulness | 42 |
| LKM | 16 |
| mindfulness, LKM | 16 |
| chakras | 3 |
| mindfulness, LKM, vipassana | 2 |
| other, (Guidée) | 2 |
| LKM, chakras | 1 |
| LKM, chakras, transcendental | 1 |
| mindfulness, LKM, chakras | 1 |
| mindfulness, LKM, chakras, other, (amour-propre et healing) | 1 |
| mindfulness, LKM, other, (J’ai fait les 42 jours de méditation par l’étude et de la méditation seul (avec bol)) | 1 |
| mindfulness, LKM, other, (en marchant) | 1 |
| mindfulness, chakras | 1 |
| mindfulness, other, (Respiration) | 1 |
| mindfulness, transcendental | 1 |
| other, (Générique, juste se coucher et se vider l’esprit) | 1 |
| other, (Libre pensée / se-laisser-porter) | 1 |
| other, (Relaxation) | 1 |
| other, (Respiration-visualisation) | 1 |
| vipassana | 1 |
get_label(data$T3_medipractice.time) %>% cat## À quelle fréquence avez-vous médité par semaine?
data %>%
count(T3_medipractice.time, sort = TRUE)| T3_medipractice.time | n |
|---|---|
| NA | 108 |
| 10-20 min | 35 |
| < 10 min | 33 |
| 20-30 min | 10 |
| 1-2 h | 8 |
| 30-60 min | 8 |
| > 2 h | 1 |
get_label(data$T3_choice.medicomp) %>% cat## Aimeriez-vous recevoir le programme de méditation gratuitement?
data %>%
count(T3_choice.medicomp, sort = TRUE)| T3_choice.medicomp | n |
|---|---|
| Oui | 104 |
| Non | 93 |
| already received | 5 |
| NA | 1 |
Consent & follow-up.
get_label(data$T3_followup) %>% cat## Aimeriez-vous être contacté.e lorsque les résultats de l’étude seront publiés dans une revue scientifique pour pouvoir en prendre connaissance?
data %>%
count(T3_followup, sort = TRUE)| T3_followup | n |
|---|---|
| Oui | 181 |
| Non | 21 |
| NA | 1 |
get_label(data$T3_consent) %>% cat## En considérant les informations précédentes, acceptez-vous de conserver votre participation à l’étude et de nous permettre d’utiliser vos données?
data %>%
count(T3_consent, sort = TRUE)| T3_consent | n |
|---|---|
| Oui | 202 |
| NA | 1 |
Here, we explore missing data, first by item and questionnaire, then
using the visdat package, and finally using Little’s MCAR
test.
# Check for nice_na
nice_na(data, scales = c(
"T1_BSCS", "T1_BAQ", "T1_NOBAGS", "T1_attitude", "T1_dehumanization",
"T1_WHS", "T1_CLS", "T2_NOBAGS", "T2_attitude", "T2_dehumanization",
"T2_SMS5", "T2_PANAS", "T2_WHS", "T2_CLS", "T2_charity", "T3_NOBAGS",
"T3_attitude", "T3_dehumanization", "T3_WHS", "T3_CLS"))| var | items | na | cells | na_percent | na_max | na_max_percent | all_na |
|---|---|---|---|---|---|---|---|
| T1_BSCS_1:T1_BSCS_7 | 7 | 0 | 1421 | 0.00 | 0 | 0.00 | 0 |
| T1_BAQ_1:T1_BAQ_13 | 13 | 1 | 2639 | 0.04 | 1 | 7.69 | 0 |
| T1_NOBAGS.1_1:T1_NOBAGS.16_1 | 20 | 0 | 4060 | 0.00 | 0 | 0.00 | 0 |
| T1_attitude_1:T1_attitude_9 | 9 | 0 | 1827 | 0.00 | 0 | 0.00 | 0 |
| T1_dehumanization_1:T1_dehumanization_9 | 9 | 1 | 1827 | 0.05 | 1 | 11.11 | 0 |
| T1_WHS_1:T1_WHS_6 | 6 | 0 | 1218 | 0.00 | 0 | 0.00 | 0 |
| T1_CLS_1:T1_CLS_21 | 21 | 0 | 4263 | 0.00 | 0 | 0.00 | 0 |
| T2_NOBAGS.1_1:T2_NOBAGS.16_1 | 20 | 0 | 4060 | 0.00 | 0 | 0.00 | 0 |
| T2_attitude_1:T2_attitude_9 | 9 | 0 | 1827 | 0.00 | 0 | 0.00 | 0 |
| T2_dehumanization_1:T2_dehumanization_9 | 9 | 1 | 1827 | 0.05 | 1 | 11.11 | 0 |
| T2_SMS5_1:T2_SMS5_6 | 6 | 0 | 1218 | 0.00 | 0 | 0.00 | 0 |
| T2_PANAS_1:T2_PANAS_10 | 10 | 0 | 2030 | 0.00 | 0 | 0.00 | 0 |
| T2_WHS_1:T2_WHS_6 | 6 | 0 | 1218 | 0.00 | 0 | 0.00 | 0 |
| T2_CLS_1:T2_CLS_21 | 21 | 0 | 4263 | 0.00 | 0 | 0.00 | 0 |
| T2_charity.moisson1_1:T2_charity.armee2_1 | 48 | 62 | 9744 | 0.64 | 22 | 45.83 | 0 |
| T3_NOBAGS.1_1:T3_NOBAGS.16_1 | 20 | 0 | 4060 | 0.00 | 0 | 0.00 | 0 |
| T3_attitude_1:T3_attitude_9 | 9 | 0 | 1827 | 0.00 | 0 | 0.00 | 0 |
| T3_dehumanization_1:T3_dehumanization_9 | 9 | 0 | 1827 | 0.00 | 0 | 0.00 | 0 |
| T3_WHS_1:T3_WHS_6 | 6 | 0 | 1218 | 0.00 | 0 | 0.00 | 0 |
| T3_CLS_1:T3_CLS_21 | 21 | 21 | 4263 | 0.49 | 21 | 100.00 | 1 |
| Total | 389 | 6763 | 78967 | 8.56 | 89 | 22.88 | 0 |
A few items are missing here and there.
Let’s check for patterns of missing data.
# Smaller subset of data for easier inspection
data %>%
# select(manualworkerId:att_check2_raw,
# condition:condition_dum) %>%
vis_miss# Let's use Little's MCAR test to confirm
# We have to proceed by "scale" because the function can only
# support 30 variables max at a time
library(naniar)## Warning: package 'naniar' was built under R version 4.2.3
# We only check for the variable that had missing data, charity
# Have to divide this one in two because it is too large for the function
data %>%
select(T2_charity.moisson1_1:T2_charity.suzuki2_1) %>%
mcar_test| statistic | df | p.value | missing.patterns |
|---|---|---|---|
| 214.1051 | 181 | 0.0466264 | 8 |
data %>%
select(T2_charity.conserv1_1:T2_charity.armee2_1) %>%
mcar_test| statistic | df | p.value | missing.patterns |
|---|---|---|---|
| 69.62351 | 61 | 0.2099961 | 6 |
Here, we impute missing data with the missForest
package, as it is one of the best imputation methods.
Why impute the data? van Ginkel explains,
Regardless of the missingness mechanism, multiple imputation is always to be preferred over listwise deletion. Under MCAR it is preferred because it results in more statistical power, under MAR it is preferred because besides more power it will give unbiased results whereas listwise deletion may not, and under NMAR it is also the preferred method because it will give less biased results than listwise deletion.
van Ginkel, J. R., Linting, M., Rippe, R. C. A., & van der Voort, A. (2020). Rebutting existing misconceptions about multiple imputation as a method for handling missing data. Journal of Personality Assessment, 102(3), 297-308. https://doi.org/10.1080/00223891.2018.1530680
Why missForest? It outperforms other imputation methods,
including the popular MICE (multiple imputation by chained equations).
You also don’t end up with several datasets, which makes it easier for
following analyses. Finally, it can be applied to mixed data types
(missings in numeric & categorical variables).
Waljee, A. K., Mukherjee, A., Singal, A. G., Zhang, Y., Warren, J., Balis, U., … & Higgins, P. D. (2013). Comparison of imputation methods for missing laboratory data in medicine. BMJ open, 3(8), e002847. https://doi.org/10.1093/bioinformatics/btr597
Stekhoven, D. J., & Bühlmann, P. (2012). MissForest—non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1), 112-118. https://doi.org/10.1093/bioinformatics/btr597
Now that we have imputed the missing data, we are ready to calculate our scale means. After reversing our items, we can then get the alphas and omegas for our different scales.
# Reverse code items 2, 4, 6, 7
data <- data %>%
mutate(across(contains("BSCS"), .names = "{col}r"))
data <- data %>%
mutate(across(ends_with(paste("BSCS", c(2, 4, 6, 7), sep = "_")),
~nice_reverse(.x, 5), .names = "{col}r"))
# Get mean BSCS
data <- data %>%
mutate(T1_BSCS = rowMeans(pick(T1_BSCS_1r:T1_BSCS_7r)))
# Get alpha & omega
data %>%
select(T1_BSCS_1r:T1_BSCS_7r) %>%
omega(nfactors = 1)## Loading required namespace: GPArotation
## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.75
## G.6: 0.75
## Omega Hierarchical: 0.75
## Omega H asymptotic: 1
## Omega Total 0.75
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_BSCS_1r 0.54 0.29 0.71 1
## T1_BSCS_2r 0.55 0.30 0.70 1
## T1_BSCS_3r 0.69 0.47 0.53 1
## T1_BSCS_4r 0.59 0.35 0.65 1
## T1_BSCS_5r 0.48 0.23 0.77 1
## T1_BSCS_6r 0.61 0.37 0.63 1
## T1_BSCS_7r 0.37 0.14 0.86 1
##
## With Sums of squares of:
## g F1*
## 2.1 0.0
##
## general/max 25873775754603412 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 14 and the fit is 0.39
## The number of observations was 203 with Chi Square = 77.58 with prob < 0.000000000079
## The root mean square of the residuals is 0.09
## The df corrected root mean square of the residuals is 0.12
## RMSEA index = 0.149 and the 10 % confidence intervals are 0.118 0.183
## BIC = 3.19
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 14 and the fit is 0.39
## The number of observations was 203 with Chi Square = 77.58 with prob < 0.000000000079
## The root mean square of the residuals is 0.09
## The df corrected root mean square of the residuals is 0.12
##
## RMSEA index = 0.149 and the 10 % confidence intervals are 0.118 0.183
## BIC = 3.19
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.88 0
## Multiple R square of scores with factors 0.77 0
## Minimum correlation of factor score estimates 0.54 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.75 0.75
## Omega general for total scores and subscales 0.75 0.75
## Omega group for total scores and subscales 0.00 0.00
# Reverse code item 7
data <- data %>%
mutate(across(contains("BAQ"), .names = "{col}r"))
data <- data %>%
mutate(across(T1_BAQ_7, ~nice_reverse(.x, 7), .names = "{col}r"))
# Get mean BAQ
data <- data %>%
mutate(T1_BAQ = rowMeans(pick(T1_BAQ_1r:T1_BAQ_12r)))
# Get alpha & omega
data %>%
select(T1_BAQ_1r:T1_BAQ_12r) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.78
## G.6: 0.82
## Omega Hierarchical: 0.79
## Omega H asymptotic: 1
## Omega Total 0.79
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_BAQ_1r 0.65 0.42 0.58 1
## T1_BAQ_2r 0.57 0.32 0.68 1
## T1_BAQ_3r 0.56 0.31 0.69 1
## T1_BAQ_4r 0.03 0.97 1
## T1_BAQ_5r 0.56 0.31 0.69 1
## T1_BAQ_6r 0.46 0.21 0.79 1
## T1_BAQ_7r 0.41 0.17 0.83 1
## T1_BAQ_8r 0.62 0.39 0.61 1
## T1_BAQ_9r 0.72 0.52 0.48 1
## T1_BAQ_10r 0.33 0.11 0.89 1
## T1_BAQ_11r 0.45 0.20 0.80 1
## T1_BAQ_12r 0.28 0.08 0.92 1
##
## With Sums of squares of:
## g F1*
## 3.1 0.0
##
## general/max 36926570500479704 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 54 and the fit is 1.38
## The number of observations was 203 with Chi Square = 270.72 with prob < 0.0000000000000000000000000000013
## The root mean square of the residuals is 0.12
## The df corrected root mean square of the residuals is 0.13
## RMSEA index = 0.141 and the 10 % confidence intervals are 0.125 0.158
## BIC = -16.19
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 54 and the fit is 1.38
## The number of observations was 203 with Chi Square = 270.72 with prob < 0.0000000000000000000000000000013
## The root mean square of the residuals is 0.12
## The df corrected root mean square of the residuals is 0.13
##
## RMSEA index = 0.141 and the 10 % confidence intervals are 0.125 0.158
## BIC = -16.19
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.91 0
## Multiple R square of scores with factors 0.83 0
## Minimum correlation of factor score estimates 0.66 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.79 0.79
## Omega general for total scores and subscales 0.79 0.79
## Omega group for total scores and subscales 0.00 0.00
data <- data %>%
mutate(T1_attitude = rowMeans(pick(T1_attitude_1:T1_attitude_9)),
T2_attitude = rowMeans(pick(T2_attitude_1:T2_attitude_9)),
T3_attitude = rowMeans(pick(T3_attitude_1:T3_attitude_9)))
# Get alpha & omega
data %>%
select(T1_attitude_1:T1_attitude_9) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.93
## G.6: 0.94
## Omega Hierarchical: 0.94
## Omega H asymptotic: 1
## Omega Total 0.94
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_attitude_1 0.92 0.84 0.16 1
## T1_attitude_2 0.78 0.61 0.39 1
## T1_attitude_3 0.91 0.82 0.18 1
## T1_attitude_4 0.85 0.73 0.27 1
## T1_attitude_5 0.83 0.68 0.32 1
## T1_attitude_6 0.83 0.68 0.32 1
## T1_attitude_7 0.39 0.15 0.85 1
## T1_attitude_8 0.83 0.69 0.31 1
## T1_attitude_9 0.69 0.47 0.53 1
##
## With Sums of squares of:
## g F1*
## 5.7 0.0
##
## general/max Inf max/min = NaN
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 27 and the fit is 0.73
## The number of observations was 203 with Chi Square = 143.81 with prob < 0.0000000000000000068
## The root mean square of the residuals is 0.05
## The df corrected root mean square of the residuals is 0.06
## RMSEA index = 0.146 and the 10 % confidence intervals are 0.123 0.17
## BIC = 0.35
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 0.73
## The number of observations was 203 with Chi Square = 143.81 with prob < 0.0000000000000000068
## The root mean square of the residuals is 0.05
## The df corrected root mean square of the residuals is 0.06
##
## RMSEA index = 0.146 and the 10 % confidence intervals are 0.123 0.17
## BIC = 0.35
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.98 0
## Multiple R square of scores with factors 0.96 0
## Minimum correlation of factor score estimates 0.91 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.94 0.94
## Omega general for total scores and subscales 0.94 0.94
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T2_attitude_1:T2_attitude_9) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## diag(.) had 0 or NA entries; non-finite result is doubtful
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.96
## G.6: 0.97
## Omega Hierarchical: 0.96
## Omega H asymptotic: 1
## Omega Total 0.96
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_attitude_1 0.92 0.85 0.15 1
## T2_attitude_2 0.87 0.75 0.25 1
## T2_attitude_3 0.94 0.89 0.11 1
## T2_attitude_4 0.88 0.77 0.23 1
## T2_attitude_5 0.91 0.82 0.18 1
## T2_attitude_6 0.91 0.83 0.17 1
## T2_attitude_7 0.56 0.32 0.68 1
## T2_attitude_8 0.93 0.86 0.14 1
## T2_attitude_9 0.84 0.71 0.29 1
##
## With Sums of squares of:
## g F1*
## 6.8 0.0
##
## general/max Inf max/min = NaN
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 27 and the fit is 0.83
## The number of observations was 203 with Chi Square = 164.1 with prob < 0.0000000000000000000013
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.03
## RMSEA index = 0.158 and the 10 % confidence intervals are 0.136 0.182
## BIC = 20.65
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 0.83
## The number of observations was 203 with Chi Square = 164.1 with prob < 0.0000000000000000000013
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.03
##
## RMSEA index = 0.158 and the 10 % confidence intervals are 0.136 0.182
## BIC = 20.65
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.99 0
## Multiple R square of scores with factors 0.98 0
## Minimum correlation of factor score estimates 0.95 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.96 0.96
## Omega general for total scores and subscales 0.96 0.96
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T3_attitude_1:T3_attitude_9) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.96
## G.6: 0.97
## Omega Hierarchical: 0.96
## Omega H asymptotic: 1
## Omega Total 0.96
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T3_attitude_1 0.96 0.93 0.07 1
## T3_attitude_2 0.83 0.68 0.32 1
## T3_attitude_3 0.94 0.88 0.12 1
## T3_attitude_4 0.86 0.74 0.26 1
## T3_attitude_5 0.92 0.85 0.15 1
## T3_attitude_6 0.88 0.78 0.22 1
## T3_attitude_7 0.61 0.37 0.63 1
## T3_attitude_8 0.90 0.81 0.19 1
## T3_attitude_9 0.79 0.62 0.38 1
##
## With Sums of squares of:
## g F1*
## 6.7 0.0
##
## general/max 119952105224445488 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 27 and the fit is 0.98
## The number of observations was 203 with Chi Square = 193.76 with prob < 0.0000000000000000000000000038
## The root mean square of the residuals is 0.04
## The df corrected root mean square of the residuals is 0.04
## RMSEA index = 0.174 and the 10 % confidence intervals are 0.152 0.198
## BIC = 50.31
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 0.98
## The number of observations was 203 with Chi Square = 193.76 with prob < 0.0000000000000000000000000038
## The root mean square of the residuals is 0.04
## The df corrected root mean square of the residuals is 0.04
##
## RMSEA index = 0.174 and the 10 % confidence intervals are 0.152 0.198
## BIC = 50.31
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.99 0
## Multiple R square of scores with factors 0.98 0
## Minimum correlation of factor score estimates 0.96 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.96 0.96
## Omega general for total scores and subscales 0.96 0.96
## Omega group for total scores and subscales 0.00 0.00
data <- data %>%
mutate(T1_dehumanization = rowMeans(pick(T1_dehumanization_1:T1_dehumanization_9)),
T2_dehumanization = rowMeans(pick(T2_dehumanization_1:T2_dehumanization_9)),
T3_dehumanization = rowMeans(pick(T3_dehumanization_1:T3_dehumanization_9)))
# Get alpha & omega
data %>%
select(T1_dehumanization_1:T1_dehumanization_9) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.94
## G.6: 0.96
## Omega Hierarchical: 0.96
## Omega H asymptotic: 1
## Omega Total 0.96
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_dehumanization_1 0.91 0.83 0.17 1
## T1_dehumanization_2 0.93 0.87 0.13 1
## T1_dehumanization_3 0.93 0.86 0.14 1
## T1_dehumanization_4 0.93 0.86 0.14 1
## T1_dehumanization_5 0.95 0.90 0.10 1
## T1_dehumanization_6 0.87 0.76 0.24 1
## T1_dehumanization_7 0.00 1.00 1
## T1_dehumanization_8 0.92 0.85 0.15 1
## T1_dehumanization_9 0.81 0.66 0.34 1
##
## With Sums of squares of:
## g F1*
## 6.6 0.0
##
## general/max 337359029506734720 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 27 and the fit is 1.24
## The number of observations was 203 with Chi Square = 244.36 with prob < 0.00000000000000000000000000000000000069
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.04
## RMSEA index = 0.199 and the 10 % confidence intervals are 0.177 0.223
## BIC = 100.9
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 1.24
## The number of observations was 203 with Chi Square = 244.36 with prob < 0.00000000000000000000000000000000000069
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.04
##
## RMSEA index = 0.199 and the 10 % confidence intervals are 0.177 0.223
## BIC = 100.9
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.99 0
## Multiple R square of scores with factors 0.98 0
## Minimum correlation of factor score estimates 0.96 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.96 0.96
## Omega general for total scores and subscales 0.96 0.96
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T2_dehumanization_1:T2_dehumanization_9) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.96
## G.6: 0.97
## Omega Hierarchical: 0.97
## Omega H asymptotic: 1
## Omega Total 0.97
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_dehumanization_1 0.98 0.96 0.04 1
## T2_dehumanization_2 0.94 0.89 0.11 1
## T2_dehumanization_3 0.94 0.89 0.11 1
## T2_dehumanization_4 0.92 0.84 0.16 1
## T2_dehumanization_5 0.96 0.92 0.08 1
## T2_dehumanization_6 0.92 0.85 0.15 1
## T2_dehumanization_7 0.00 1.00 1
## T2_dehumanization_8 0.95 0.91 0.09 1
## T2_dehumanization_9 0.94 0.88 0.12 1
##
## With Sums of squares of:
## g F1*
## 7.1 0.0
##
## general/max 391733323604288832 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 27 and the fit is 1.09
## The number of observations was 203 with Chi Square = 215.86 with prob < 0.00000000000000000000000000000023
## The root mean square of the residuals is 0.02
## The df corrected root mean square of the residuals is 0.02
## RMSEA index = 0.186 and the 10 % confidence intervals are 0.163 0.21
## BIC = 72.41
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 1.09
## The number of observations was 203 with Chi Square = 215.86 with prob < 0.00000000000000000000000000000023
## The root mean square of the residuals is 0.02
## The df corrected root mean square of the residuals is 0.02
##
## RMSEA index = 0.186 and the 10 % confidence intervals are 0.163 0.21
## BIC = 72.41
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.99 0
## Multiple R square of scores with factors 0.99 0
## Minimum correlation of factor score estimates 0.98 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.97 0.97
## Omega general for total scores and subscales 0.97 0.97
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T3_dehumanization_1:T3_dehumanization_9) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.95
## G.6: 0.97
## Omega Hierarchical: 0.96
## Omega H asymptotic: 1
## Omega Total 0.96
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T3_dehumanization_1 0.95 0.90 0.10 1
## T3_dehumanization_2 0.93 0.87 0.13 1
## T3_dehumanization_3 0.93 0.87 0.13 1
## T3_dehumanization_4 0.95 0.91 0.09 1
## T3_dehumanization_5 0.95 0.91 0.09 1
## T3_dehumanization_6 0.89 0.78 0.22 1
## T3_dehumanization_7 0.01 0.99 1
## T3_dehumanization_8 0.93 0.86 0.14 1
## T3_dehumanization_9 0.84 0.70 0.30 1
##
## With Sums of squares of:
## g F1*
## 6.8 0.0
##
## general/max Inf max/min = NaN
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 27 and the fit is 2.32
## The number of observations was 203 with Chi Square = 457.5 with prob < 0.000000000000000000000000000000000000000000000000000000000000000000000000000000087
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.04
## RMSEA index = 0.28 and the 10 % confidence intervals are 0.259 0.304
## BIC = 314.04
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 2.32
## The number of observations was 203 with Chi Square = 457.5 with prob < 0.000000000000000000000000000000000000000000000000000000000000000000000000000000087
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.04
##
## RMSEA index = 0.28 and the 10 % confidence intervals are 0.259 0.304
## BIC = 314.04
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.99 0
## Multiple R square of scores with factors 0.98 0
## Minimum correlation of factor score estimates 0.96 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.96 0.96
## Omega general for total scores and subscales 0.96 0.96
## Omega group for total scores and subscales 0.00 0.00
data <- data %>%
mutate(across(contains("NOBAGS"), .names = "{col}r"))
# Reverse code NOBAGS (items 1:2, 5:6, 10,12, 14:16, 20)
data <- data %>%
mutate(across(ends_with(paste0("NOBAGS.", c(
"1_1", "1_2", "3_1", "3_2", "6_1", "8_1", "10_1", "12_1", "16_1"))),
~nice_reverse(.x, 4), .names = "{col}r"))
# Get mean NOBAGS
data <- data %>%
mutate(T1_NOBAGS = rowMeans(pick(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r)),
T2_NOBAGS = rowMeans(pick(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r)),
T3_NOBAGS = rowMeans(pick(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r)))
# Get alpha & omega
data %>%
select(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.85
## G.6: 0.91
## Omega Hierarchical: 0.84
## Omega H asymptotic: 0.99
## Omega Total 0.85
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_NOBAGS.1_1r 0.63 0.40 0.60 1
## T1_NOBAGS.1_2r 0.52 0.27 0.73 1
## T1_NOBAGS.2_1r 0.53 0.28 0.72 1
## T1_NOBAGS.2_2r 0.46 0.21 0.79 1
## T1_NOBAGS.3_1r 0.64 0.41 0.59 1
## T1_NOBAGS.3_2r 0.50 0.25 0.75 1
## T1_NOBAGS.4_1r 0.67 0.45 0.55 1
## T1_NOBAGS.4_2r 0.47 0.22 0.78 1
## T1_NOBAGS.5_1r 0.51 0.26 0.74 1
## T1_NOBAGS.6_1r 0.42 0.18 0.82 1
## T1_NOBAGS.7_1r 0.57 0.33 0.67 1
## T1_NOBAGS.8_1r 0.45 0.20 0.80 1
## T1_NOBAGS.9_1r 0.44 0.20 0.80 1
## T1_NOBAGS.10_1r 0.34 0.12 0.88 1
## T1_NOBAGS.11_1r- 0.27 0.07 0.93 1
## T1_NOBAGS.12_1r 0.30 0.09 0.91 1
## T1_NOBAGS.13_1r 0.43 0.18 0.82 1
## T1_NOBAGS.14_1r 0.45 0.21 0.79 1
## T1_NOBAGS.15_1r 0.47 0.22 0.78 1
## T1_NOBAGS.16_1r 0.03 0.97 1
##
## With Sums of squares of:
## g F1*
## 4.6 0.0
##
## general/max 14975690996060718 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 170 and the fit is 6.58
## The number of observations was 203 with Chi Square = 1276.17 with prob < 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is 0.16
## The df corrected root mean square of the residuals is 0.17
## RMSEA index = 0.179 and the 10 % confidence intervals are 0.17 0.189
## BIC = 372.93
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 170 and the fit is 6.58
## The number of observations was 203 with Chi Square = 1276.17 with prob < 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is 0.16
## The df corrected root mean square of the residuals is 0.17
##
## RMSEA index = 0.179 and the 10 % confidence intervals are 0.17 0.189
## BIC = 372.93
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.93 0
## Multiple R square of scores with factors 0.87 0
## Minimum correlation of factor score estimates 0.73 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.85 0.84
## Omega general for total scores and subscales 0.84 0.84
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.88
## G.6: 0.94
## Omega Hierarchical: 0.88
## Omega H asymptotic: 0.99
## Omega Total 0.88
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_NOBAGS.1_1r 0.53 0.28 0.72 1
## T2_NOBAGS.1_2r 0.41 0.17 0.83 1
## T2_NOBAGS.2_1r 0.68 0.46 0.54 1
## T2_NOBAGS.2_2r 0.69 0.47 0.53 1
## T2_NOBAGS.3_1r 0.63 0.40 0.60 1
## T2_NOBAGS.3_2r 0.54 0.29 0.71 1
## T2_NOBAGS.4_1r 0.64 0.41 0.59 1
## T2_NOBAGS.4_2r 0.61 0.37 0.63 1
## T2_NOBAGS.5_1r 0.68 0.47 0.53 1
## T2_NOBAGS.6_1r 0.45 0.20 0.80 1
## T2_NOBAGS.7_1r 0.68 0.46 0.54 1
## T2_NOBAGS.8_1r 0.42 0.18 0.82 1
## T2_NOBAGS.9_1r 0.46 0.21 0.79 1
## T2_NOBAGS.10_1r 0.37 0.13 0.87 1
## T2_NOBAGS.11_1r- 0.38 0.15 0.85 1
## T2_NOBAGS.12_1r 0.46 0.21 0.79 1
## T2_NOBAGS.13_1r 0.47 0.22 0.78 1
## T2_NOBAGS.14_1r 0.47 0.22 0.78 1
## T2_NOBAGS.15_1r 0.42 0.18 0.82 1
## T2_NOBAGS.16_1r 0.42 0.18 0.82 1
##
## With Sums of squares of:
## g F1*
## 5.7 0.0
##
## general/max 17017814006844390 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 170 and the fit is 7.91
## The number of observations was 203 with Chi Square = 1532.35 with prob < 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012
## The root mean square of the residuals is 0.18
## The df corrected root mean square of the residuals is 0.19
## RMSEA index = 0.199 and the 10 % confidence intervals are 0.19 0.208
## BIC = 629.1
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 170 and the fit is 7.91
## The number of observations was 203 with Chi Square = 1532.35 with prob < 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012
## The root mean square of the residuals is 0.18
## The df corrected root mean square of the residuals is 0.19
##
## RMSEA index = 0.199 and the 10 % confidence intervals are 0.19 0.208
## BIC = 629.1
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.95 0
## Multiple R square of scores with factors 0.90 0
## Minimum correlation of factor score estimates 0.79 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.88 0.88
## Omega general for total scores and subscales 0.88 0.88
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.89
## G.6: 0.94
## Omega Hierarchical: 0.89
## Omega H asymptotic: 0.99
## Omega Total 0.89
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T3_NOBAGS.1_1r 0.63 0.40 0.60 1
## T3_NOBAGS.1_2r 0.44 0.20 0.80 1
## T3_NOBAGS.2_1r 0.75 0.56 0.44 1
## T3_NOBAGS.2_2r 0.69 0.48 0.52 1
## T3_NOBAGS.3_1r 0.64 0.41 0.59 1
## T3_NOBAGS.3_2r 0.55 0.30 0.70 1
## T3_NOBAGS.4_1r 0.72 0.52 0.48 1
## T3_NOBAGS.4_2r 0.61 0.37 0.63 1
## T3_NOBAGS.5_1r 0.68 0.46 0.54 1
## T3_NOBAGS.6_1r 0.48 0.23 0.77 1
## T3_NOBAGS.7_1r 0.68 0.46 0.54 1
## T3_NOBAGS.8_1r 0.43 0.19 0.81 1
## T3_NOBAGS.9_1r 0.51 0.26 0.74 1
## T3_NOBAGS.10_1r 0.42 0.18 0.82 1
## T3_NOBAGS.11_1r- 0.35 0.12 0.88 1
## T3_NOBAGS.12_1r 0.42 0.18 0.82 1
## T3_NOBAGS.13_1r 0.46 0.21 0.79 1
## T3_NOBAGS.14_1r 0.54 0.29 0.71 1
## T3_NOBAGS.15_1r 0.47 0.22 0.78 1
## T3_NOBAGS.16_1r 0.34 0.12 0.88 1
##
## With Sums of squares of:
## g F1*
## 6.1 0.0
##
## general/max 14727181914176212 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 170 and the fit is 7.67
## The number of observations was 203 with Chi Square = 1487.54 with prob < 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000052
## The root mean square of the residuals is 0.17
## The df corrected root mean square of the residuals is 0.18
## RMSEA index = 0.195 and the 10 % confidence intervals are 0.187 0.205
## BIC = 584.29
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 170 and the fit is 7.67
## The number of observations was 203 with Chi Square = 1487.54 with prob < 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000052
## The root mean square of the residuals is 0.17
## The df corrected root mean square of the residuals is 0.18
##
## RMSEA index = 0.195 and the 10 % confidence intervals are 0.187 0.205
## BIC = 584.29
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.95 0
## Multiple R square of scores with factors 0.91 0
## Minimum correlation of factor score estimates 0.82 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.89 0.89
## Omega general for total scores and subscales 0.89 0.89
## Omega group for total scores and subscales 0.00 0.00
data <- data %>%
mutate(T1_WHS = rowMeans(pick(T1_WHS_1:T1_WHS_6)),
T2_WHS = rowMeans(pick(T2_WHS_1:T2_WHS_6)),
T3_WHS = rowMeans(pick(T3_WHS_1:T3_WHS_6)))
# Get alpha & omega
data %>%
select(T1_WHS_1:T1_WHS_6) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.65
## G.6: 0.66
## Omega Hierarchical: 0.64
## Omega H asymptotic: 0.96
## Omega Total 0.67
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_WHS_1 0.72 0.52 0.48 1
## T1_WHS_2 0.57 0.32 0.68 1
## T1_WHS_3 0.24 0.06 0.94 1
## T1_WHS_4 0.46 0.21 0.79 1
## T1_WHS_5 0.25 0.06 0.94 1
## T1_WHS_6 0.67 0.45 0.55 1
##
## With Sums of squares of:
## g F1*
## 1.6 0.0
##
## general/max 14564732675507748 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 9 and the fit is 0.26
## The number of observations was 203 with Chi Square = 50.63 with prob < 0.000000082
## The root mean square of the residuals is 0.11
## The df corrected root mean square of the residuals is 0.15
## RMSEA index = 0.151 and the 10 % confidence intervals are 0.112 0.193
## BIC = 2.82
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9 and the fit is 0.26
## The number of observations was 203 with Chi Square = 50.63 with prob < 0.000000082
## The root mean square of the residuals is 0.11
## The df corrected root mean square of the residuals is 0.15
##
## RMSEA index = 0.151 and the 10 % confidence intervals are 0.112 0.193
## BIC = 2.82
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.86 0
## Multiple R square of scores with factors 0.73 0
## Minimum correlation of factor score estimates 0.46 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.67 0.64
## Omega general for total scores and subscales 0.64 0.64
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T2_WHS_1:T2_WHS_6) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.67
## G.6: 0.67
## Omega Hierarchical: 0.67
## Omega H asymptotic: 0.98
## Omega Total 0.68
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_WHS_1 0.68 0.47 0.53 1
## T2_WHS_2 0.60 0.36 0.64 1
## T2_WHS_3 0.38 0.14 0.86 1
## T2_WHS_4 0.41 0.16 0.84 1
## T2_WHS_5 0.34 0.11 0.89 1
## T2_WHS_6 0.62 0.39 0.61 1
##
## With Sums of squares of:
## g F1*
## 1.6 0.0
##
## general/max 58885651346390504 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 9 and the fit is 0.26
## The number of observations was 203 with Chi Square = 51.27 with prob < 0.000000062
## The root mean square of the residuals is 0.11
## The df corrected root mean square of the residuals is 0.14
## RMSEA index = 0.152 and the 10 % confidence intervals are 0.113 0.194
## BIC = 3.45
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9 and the fit is 0.26
## The number of observations was 203 with Chi Square = 51.27 with prob < 0.000000062
## The root mean square of the residuals is 0.11
## The df corrected root mean square of the residuals is 0.14
##
## RMSEA index = 0.152 and the 10 % confidence intervals are 0.113 0.194
## BIC = 3.45
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.85 0
## Multiple R square of scores with factors 0.72 0
## Minimum correlation of factor score estimates 0.43 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.68 0.67
## Omega general for total scores and subscales 0.67 0.67
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T3_WHS_1:T3_WHS_6) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.69
## G.6: 0.69
## Omega Hierarchical: 0.69
## Omega H asymptotic: 0.99
## Omega Total 0.69
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T3_WHS_1 0.50 0.25 0.75 1
## T3_WHS_2 0.57 0.33 0.67 1
## T3_WHS_3 0.37 0.14 0.86 1
## T3_WHS_4 0.45 0.20 0.80 1
## T3_WHS_5 0.50 0.25 0.75 1
## T3_WHS_6 0.71 0.50 0.50 1
##
## With Sums of squares of:
## g F1*
## 1.7 0.0
##
## general/max 10006094449321290 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 9 and the fit is 0.33
## The number of observations was 203 with Chi Square = 65.07 with prob < 0.00000000014
## The root mean square of the residuals is 0.12
## The df corrected root mean square of the residuals is 0.16
## RMSEA index = 0.175 and the 10 % confidence intervals are 0.137 0.217
## BIC = 17.25
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9 and the fit is 0.33
## The number of observations was 203 with Chi Square = 65.07 with prob < 0.00000000014
## The root mean square of the residuals is 0.12
## The df corrected root mean square of the residuals is 0.16
##
## RMSEA index = 0.175 and the 10 % confidence intervals are 0.137 0.217
## BIC = 17.25
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.85 0
## Multiple R square of scores with factors 0.72 0
## Minimum correlation of factor score estimates 0.44 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.69 0.69
## Omega general for total scores and subscales 0.69 0.69
## Omega group for total scores and subscales 0.00 0.00
data <- data %>%
mutate(T1_CLS = rowMeans(pick(T1_CLS_1:T1_CLS_21)),
T2_CLS = rowMeans(pick(T2_CLS_1:T2_CLS_21)),
T3_CLS = rowMeans(pick(T3_CLS_1:T3_CLS_21)))
# Get alpha & omega
data %>%
select(T1_CLS_1:T1_CLS_21) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.94
## G.6: 0.95
## Omega Hierarchical: 0.94
## Omega H asymptotic: 1
## Omega Total 0.94
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T1_CLS_1 0.67 0.44 0.56 1
## T1_CLS_2 0.64 0.41 0.59 1
## T1_CLS_3 0.74 0.55 0.45 1
## T1_CLS_4 0.65 0.42 0.58 1
## T1_CLS_5 0.71 0.50 0.50 1
## T1_CLS_6 0.74 0.55 0.45 1
## T1_CLS_7 0.54 0.29 0.71 1
## T1_CLS_8 0.72 0.51 0.49 1
## T1_CLS_9 0.77 0.59 0.41 1
## T1_CLS_10 0.67 0.45 0.55 1
## T1_CLS_11 0.60 0.35 0.65 1
## T1_CLS_12 0.72 0.52 0.48 1
## T1_CLS_13 0.42 0.18 0.82 1
## T1_CLS_14 0.53 0.28 0.72 1
## T1_CLS_15 0.70 0.49 0.51 1
## T1_CLS_16 0.60 0.36 0.64 1
## T1_CLS_17 0.70 0.49 0.51 1
## T1_CLS_18 0.63 0.40 0.60 1
## T1_CLS_19 0.59 0.35 0.65 1
## T1_CLS_20 0.65 0.43 0.57 1
## T1_CLS_21 0.66 0.44 0.56 1
##
## With Sums of squares of:
## g F1*
## 9 0
##
## general/max 54179614793415584 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 189 and the fit is 1.91
## The number of observations was 203 with Chi Square = 369.45 with prob < 0.000000000000089
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.06
## RMSEA index = 0.068 and the 10 % confidence intervals are 0.058 0.079
## BIC = -634.75
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 189 and the fit is 1.91
## The number of observations was 203 with Chi Square = 369.45 with prob < 0.000000000000089
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.06
##
## RMSEA index = 0.068 and the 10 % confidence intervals are 0.058 0.079
## BIC = -634.75
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.97 0
## Multiple R square of scores with factors 0.94 0
## Minimum correlation of factor score estimates 0.89 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.94 0.94
## Omega general for total scores and subscales 0.94 0.94
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T2_CLS_1:T2_CLS_21) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.94
## G.6: 0.96
## Omega Hierarchical: 0.95
## Omega H asymptotic: 1
## Omega Total 0.95
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_CLS_1 0.69 0.48 0.52 1
## T2_CLS_2 0.71 0.50 0.50 1
## T2_CLS_3 0.79 0.62 0.38 1
## T2_CLS_4 0.63 0.40 0.60 1
## T2_CLS_5 0.68 0.46 0.54 1
## T2_CLS_6 0.77 0.59 0.41 1
## T2_CLS_7 0.62 0.38 0.62 1
## T2_CLS_8 0.72 0.51 0.49 1
## T2_CLS_9 0.77 0.60 0.40 1
## T2_CLS_10 0.73 0.53 0.47 1
## T2_CLS_11 0.64 0.41 0.59 1
## T2_CLS_12 0.76 0.57 0.43 1
## T2_CLS_13 0.30 0.09 0.91 1
## T2_CLS_14 0.57 0.32 0.68 1
## T2_CLS_15 0.79 0.62 0.38 1
## T2_CLS_16 0.61 0.37 0.63 1
## T2_CLS_17 0.70 0.49 0.51 1
## T2_CLS_18 0.62 0.39 0.61 1
## T2_CLS_19 0.64 0.40 0.60 1
## T2_CLS_20 0.64 0.42 0.58 1
## T2_CLS_21 0.66 0.43 0.57 1
##
## With Sums of squares of:
## g F1*
## 9.6 0.0
##
## general/max 34494890394986312 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 189 and the fit is 2.68
## The number of observations was 203 with Chi Square = 517.68 with prob < 0.000000000000000000000000000000022
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.07
## RMSEA index = 0.092 and the 10 % confidence intervals are 0.083 0.102
## BIC = -486.51
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 189 and the fit is 2.68
## The number of observations was 203 with Chi Square = 517.68 with prob < 0.000000000000000000000000000000022
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.07
##
## RMSEA index = 0.092 and the 10 % confidence intervals are 0.083 0.102
## BIC = -486.51
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.98 0
## Multiple R square of scores with factors 0.95 0
## Minimum correlation of factor score estimates 0.90 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.95 0.95
## Omega general for total scores and subscales 0.95 0.95
## Omega group for total scores and subscales 0.00 0.00
data %>%
select(T3_CLS_1:T3_CLS_21) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.95
## G.6: 0.96
## Omega Hierarchical: 0.95
## Omega H asymptotic: 1
## Omega Total 0.95
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T3_CLS_1 0.73 0.53 0.47 1
## T3_CLS_2 0.76 0.57 0.43 1
## T3_CLS_3 0.77 0.59 0.41 1
## T3_CLS_4 0.71 0.50 0.50 1
## T3_CLS_5 0.69 0.48 0.52 1
## T3_CLS_6 0.79 0.63 0.37 1
## T3_CLS_7 0.65 0.42 0.58 1
## T3_CLS_8 0.71 0.50 0.50 1
## T3_CLS_9 0.77 0.59 0.41 1
## T3_CLS_10 0.73 0.54 0.46 1
## T3_CLS_11 0.66 0.43 0.57 1
## T3_CLS_12 0.79 0.62 0.38 1
## T3_CLS_13 0.46 0.21 0.79 1
## T3_CLS_14 0.55 0.30 0.70 1
## T3_CLS_15 0.80 0.64 0.36 1
## T3_CLS_16 0.61 0.38 0.62 1
## T3_CLS_17 0.74 0.55 0.45 1
## T3_CLS_18 0.69 0.48 0.52 1
## T3_CLS_19 0.61 0.37 0.63 1
## T3_CLS_20 0.63 0.39 0.61 1
## T3_CLS_21 0.68 0.46 0.54 1
##
## With Sums of squares of:
## g F1*
## 10 0
##
## general/max 36644823991854128 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 189 and the fit is 2.78
## The number of observations was 203 with Chi Square = 537.38 with prob < 0.000000000000000000000000000000000038
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.07
## RMSEA index = 0.095 and the 10 % confidence intervals are 0.086 0.105
## BIC = -466.81
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 189 and the fit is 2.78
## The number of observations was 203 with Chi Square = 537.38 with prob < 0.000000000000000000000000000000000038
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.07
##
## RMSEA index = 0.095 and the 10 % confidence intervals are 0.086 0.105
## BIC = -466.81
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.98 0
## Multiple R square of scores with factors 0.96 0
## Minimum correlation of factor score estimates 0.91 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.95 0.95
## Omega general for total scores and subscales 0.95 0.95
## Omega group for total scores and subscales 0.00 0.00
data <- data %>%
mutate(across(contains("SMS5"), .names = "{col}r"))
# Reverse code SMS5 (items 3 et 5)
data <- data %>%
mutate(across(ends_with(paste0("SMS5_", c(3, 5))),
~nice_reverse(.x, 5), .names = "{col}r"))
# Get mean SMS5
data <- data %>%
mutate(T2_SMS5 = rowMeans(pick(T2_SMS5_1r:T2_SMS5_6r)))
# Get alpha & omega
data %>%
select(T2_SMS5_1r:T2_SMS5_6r) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.65
## G.6: 0.65
## Omega Hierarchical: 0.67
## Omega H asymptotic: 1
## Omega Total 0.67
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_SMS5_1r 0.63 0.40 0.60 1
## T2_SMS5_2r 0.00 1.00 1
## T2_SMS5_3r 0.50 0.25 0.75 1
## T2_SMS5_4r 0.58 0.34 0.66 1
## T2_SMS5_5r 0.59 0.35 0.65 1
## T2_SMS5_6r 0.62 0.39 0.61 1
##
## With Sums of squares of:
## g F1*
## 1.7 0.0
##
## general/max 11338623069083538 max/min = 1
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 9 and the fit is 0.15
## The number of observations was 203 with Chi Square = 30.12 with prob < 0.00042
## The root mean square of the residuals is 0.07
## The df corrected root mean square of the residuals is 0.09
## RMSEA index = 0.107 and the 10 % confidence intervals are 0.067 0.151
## BIC = -17.7
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9 and the fit is 0.15
## The number of observations was 203 with Chi Square = 30.12 with prob < 0.00042
## The root mean square of the residuals is 0.07
## The df corrected root mean square of the residuals is 0.09
##
## RMSEA index = 0.107 and the 10 % confidence intervals are 0.067 0.151
## BIC = -17.7
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.85 0
## Multiple R square of scores with factors 0.73 0
## Minimum correlation of factor score estimates 0.46 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.67 0.67
## Omega general for total scores and subscales 0.67 0.67
## Omega group for total scores and subscales 0.00 0.00
# Get mean of PANAS
# Positive affect = 1, 3, 5, 7, 9
# Negative affect = 2, 4, 6, 8, 10
data <- data %>% mutate(
T2_PANAS_pos = rowMeans(pick(paste0("T2_PANAS_", seq(1, 9, 2)))),
T2_PANAS_neg = rowMeans(pick(paste0("T2_PANAS_", seq(2, 10, 2)))))
# Get alpha & omega
# Positive affect
data %>%
select(all_of(paste0("T2_PANAS_", seq(1, 9, 2)))) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.86
## G.6: 0.84
## Omega Hierarchical: 0.86
## Omega H asymptotic: 1
## Omega Total 0.86
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_PANAS_1 0.73 0.53 0.47 1
## T2_PANAS_3 0.70 0.49 0.51 1
## T2_PANAS_5 0.79 0.63 0.37 1
## T2_PANAS_7 0.67 0.45 0.55 1
## T2_PANAS_9 0.84 0.70 0.30 1
##
## With Sums of squares of:
## g F1*
## 2.8 0.0
##
## general/max Inf max/min = NaN
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 5 and the fit is 0.05
## The number of observations was 203 with Chi Square = 10.65 with prob < 0.059
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.04
## RMSEA index = 0.074 and the 10 % confidence intervals are 0 0.138
## BIC = -15.91
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 5 and the fit is 0.05
## The number of observations was 203 with Chi Square = 10.65 with prob < 0.059
## The root mean square of the residuals is 0.03
## The df corrected root mean square of the residuals is 0.04
##
## RMSEA index = 0.074 and the 10 % confidence intervals are 0 0.138
## BIC = -15.91
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.94 0
## Multiple R square of scores with factors 0.88 0
## Minimum correlation of factor score estimates 0.75 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.86 0.86
## Omega general for total scores and subscales 0.86 0.86
## Omega group for total scores and subscales 0.00 0.00
# Negative affect
data %>%
select(all_of(paste0("T2_PANAS_", seq(2, 10, 2)))) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## diag(.) had 0 or NA entries; non-finite result is doubtful
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.8
## G.6: 0.8
## Omega Hierarchical: 0.82
## Omega H asymptotic: 1
## Omega Total 0.82
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_PANAS_2 0.85 0.73 0.27 1
## T2_PANAS_4 0.34 0.11 0.89 1
## T2_PANAS_6 0.79 0.62 0.38 1
## T2_PANAS_8 0.67 0.45 0.55 1
## T2_PANAS_10 0.73 0.54 0.46 1
##
## With Sums of squares of:
## g F1*
## 2.5 0.0
##
## general/max Inf max/min = NaN
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 5 and the fit is 0.16
## The number of observations was 203 with Chi Square = 31.84 with prob < 0.0000064
## The root mean square of the residuals is 0.07
## The df corrected root mean square of the residuals is 0.1
## RMSEA index = 0.163 and the 10 % confidence intervals are 0.112 0.219
## BIC = 5.28
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 5 and the fit is 0.16
## The number of observations was 203 with Chi Square = 31.84 with prob < 0.0000064
## The root mean square of the residuals is 0.07
## The df corrected root mean square of the residuals is 0.1
##
## RMSEA index = 0.163 and the 10 % confidence intervals are 0.112 0.219
## BIC = 5.28
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.93 0
## Multiple R square of scores with factors 0.87 0
## Minimum correlation of factor score estimates 0.73 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.82 0.82
## Omega general for total scores and subscales 0.82 0.82
## Omega group for total scores and subscales 0.00 0.00
data <- data %>% mutate(
T2_Charity = rowMeans(pick(contains("charity") & ends_with("1_1"))),
T2_Familiarity = rowMeans(pick(contains("charity") & ends_with("2_1"))))
# Get alpha & omega
data %>%
select(contains("charity") & ends_with("1_1")) %>%
omega(nfactors = 1)## Omega_h for 1 factor is not meaningful, just omega_t
## Warning in schmid(m, nfactors, fm, digits, rotate = rotate, n.obs = n.obs, :
## Omega_h and Omega_asymptotic are not meaningful with one factor
## Warning in cov2cor(t(w) %*% r %*% w): diag(.) had 0 or NA entries; non-finite
## result is doubtful
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.98
## G.6: 0.99
## Omega Hierarchical: 0.98
## Omega H asymptotic: 1
## Omega Total 0.98
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* h2 u2 p2
## T2_charity.moisson1_1 0.84 0.71 0.29 1
## T2_charity.bonneau1_1 0.85 0.72 0.28 1
## T2_charity.AMDI1_1 0.86 0.74 0.26 1
## T2_charity.CAAM1_1 0.86 0.74 0.26 1
## T2_charity.centrefem1_1 0.80 0.64 0.36 1
## T2_charity.LGBTQ1_1 0.74 0.55 0.45 1
## T2_charity.equiterre1_1 0.87 0.75 0.25 1
## T2_charity.dejeuner1_1 0.80 0.65 0.35 1
## T2_charity.cancer1_1 0.87 0.75 0.25 1
## T2_charity.parkinson1_1 0.87 0.76 0.24 1
## T2_charity.OXFAM1_1 0.81 0.66 0.34 1
## T2_charity.papillon1_1 0.85 0.72 0.28 1
## T2_charity.croix1_1 0.86 0.74 0.26 1
## T2_charity.centraide1_1 0.87 0.76 0.24 1
## T2_charity.autisme1_1 0.82 0.68 0.32 1
## T2_charity.suzuki1_1 0.78 0.62 0.38 1
## T2_charity.conserv1_1 0.81 0.65 0.35 1
## T2_charity.coeur1_1 0.86 0.73 0.27 1
## T2_charity.UNICEF1_1 0.82 0.68 0.32 1
## T2_charity.amnistie1_1 0.81 0.66 0.34 1
## T2_charity.green1_1 0.77 0.59 0.41 1
## T2_charity.WWF1_1 0.79 0.62 0.38 1
## T2_charity.MSF1_1 0.86 0.74 0.26 1
## T2_charity.armee1_1 0.75 0.57 0.43 1
##
## With Sums of squares of:
## g F1*
## 16 0
##
## general/max Inf max/min = NaN
## mean percent general = 1 with sd = 0 and cv of 0
## Explained Common Variance of the general factor = 1
##
## The degrees of freedom are 252 and the fit is 4.94
## The number of observations was 203 with Chi Square = 950.72 with prob < 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is 0.05
## The df corrected root mean square of the residuals is 0.05
## RMSEA index = 0.117 and the 10 % confidence intervals are 0.109 0.125
## BIC = -388.21
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 252 and the fit is 4.94
## The number of observations was 203 with Chi Square = 950.72 with prob < 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000011
## The root mean square of the residuals is 0.05
## The df corrected root mean square of the residuals is 0.05
##
## RMSEA index = 0.117 and the 10 % confidence intervals are 0.109 0.125
## BIC = -388.21
##
## Measures of factor score adequacy
## g F1*
## Correlation of scores with factors 0.99 0
## Multiple R square of scores with factors 0.98 0
## Minimum correlation of factor score estimates 0.96 -1
##
## Total, General and Subset omega for each subset
## g F1*
## Omega total for total scores and subscales 0.98 0.98
## Omega general for total scores and subscales 0.98 0.98
## Omega group for total scores and subscales 0.00 0.00
In this section, we: (a) test assumptions of normality, (b) transform variables violating assumptions, (c) test assumptions of homoscedasticity, and (d) identify and winsorize outliers.
At this stage, we define a list of our relevant variables.
# Make list of DVs
col.list <- c("T1_blastintensity", "T1_blastduration", "T1_blastintensity.duration",
"T2_blastintensity", "T2_blastduration", "T2_blastintensity.duration",
"T3_blastintensity", "T3_blastduration", "T3_blastintensity.duration",
"T1_BSCS", "T1_BAQ", "T1_attitude", "T2_attitude", "T3_attitude",
"T1_dehumanization", "T2_dehumanization", "T3_dehumanization",
"T1_NOBAGS", "T1_WHS", "T2_WHS", "T3_WHS", "T1_CLS", "T2_CLS", "T3_CLS",
"T2_SMS5", "T2_PANAS_pos", "T2_PANAS_neg", "T2_Charity",
"T2_Familiarity", "T2_memory.altruistic", "T2_memory.aggressive")lapply(col.list, function(x)
nice_normality(data,
variable = x,
title = x,
group = "T1_Group",
shapiro = TRUE,
histogram = TRUE))## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
##
## [[15]]
##
## [[16]]
##
## [[17]]
##
## [[18]]
##
## [[19]]
##
## [[20]]
##
## [[21]]
##
## [[22]]
##
## [[23]]
##
## [[24]]
##
## [[25]]
##
## [[26]]
##
## [[27]]
##
## [[28]]
##
## [[29]]
##
## [[30]]
##
## [[31]]
Several variables are clearly skewed. Let’s apply transformations. But first, let’s deal with the working memory task, SOPT (Self-Ordered Pointing Task). It is clearly problematic.
The function below transforms variables according to the best
possible transformation (via the bestNormalize package),
and also standardizes the variables.
Note. The I(x) transformations above are actually not transformations, but a shorthand function for passing the data “as is”. Suggesting the package estimated the various attempted transformations did not improve normality in those cases, so no transformation is used. This only appears when standardize is set to FALSE. When set to TRUE, for those variables, it is actually center_scale(x), suggesting that the data are only CENTERED because they need no transformation (no need to be scaled), only to be centered.
Let’s check if normality was corrected.
Looks rather reasonable now, though not perfect (fortunately t-tests are quite robust against violations of normality).
We can now resume with the next step: checking variance.
# Plotting variance
plots(lapply(col.list, function(x) {
nice_varplot(data, x, group = "T1_Group")
}),
n_columns = 2)Variance looks good. No group has four times the variance of any other group. We can now resume with checking outliers.
We check outliers visually with the plot_outliers
function, which draws red lines at +/- 3 median absolute deviations.
plots(lapply(col.list, function(x) {
plot_outliers(data, x, group = "T1_Group", ytitle = x, binwidth = 0.3)
}),
n_columns = 2)There are some outliers, but nothing unreasonable. Let’s still check with the 3 median absolute deviations (MAD) method.
data %>%
as.data.frame %>%
filter(T1_Group == "Waitlist") %>%
find_mad(col.list, criteria = 3)## 39 outlier(s) based on 3 median absolute deviations for variable(s):
## T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T1_dehumanization, T2_dehumanization, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive
##
## The following participants were considered outliers for more than one variable:
##
## Row n
## 1 6 6
## 2 7 5
## 3 16 2
## 4 20 3
## 5 26 2
## 6 37 3
## 7 43 2
## 8 49 2
## 9 50 2
## 10 51 3
## 11 66 2
## 12 69 3
## 13 70 2
## 14 83 2
## 15 87 4
##
## Outliers per variable:
##
## $T1_blastintensity.duration
## Row T1_blastintensity.duration_mad
## 1 37 4.279092
##
## $T3_blastintensity.duration
## Row T3_blastintensity.duration_mad
## 1 4 3.432722
## 2 9 3.353952
##
## $T1_attitude
## Row T1_attitude_mad
## 1 6 -5.178348
## 2 15 -3.054790
## 3 26 -3.220149
## 4 50 -3.098306
##
## $T2_attitude
## Row T2_attitude_mad
## 1 6 -5.750247
## 2 26 -3.052284
## 3 66 -3.052284
##
## $T3_attitude
## Row T3_attitude_mad
## 1 6 -6.397443
##
## $T1_dehumanization
## Row T1_dehumanization_mad
## 1 5 -8.408651
## 2 6 -13.040155
## 3 7 -5.395926
## 4 16 -7.329466
## 5 20 -6.542560
## 6 41 -3.934529
## 7 43 -6.902289
## 8 50 -3.709699
## 9 51 -5.440892
## 10 58 -5.845587
## 11 64 -3.507352
## 12 69 -6.025451
## 13 79 -3.147624
## 14 83 -4.721435
## 15 87 -5.575790
## 16 88 -4.609020
##
## $T2_dehumanization
## Row T2_dehumanization_mad
## 1 6 -13.799716
## 2 7 -3.919338
## 3 20 -5.541762
## 4 37 -3.208388
## 5 51 -3.500060
## 6 69 -3.554749
## 7 70 -4.265698
## 8 87 -5.049566
##
## $T3_dehumanization
## Row T3_dehumanization_mad
## 1 6 -13.267943
## 2 7 -3.594326
## 3 16 -4.393065
## 4 20 -7.304025
## 5 37 -3.008584
## 6 43 -3.416828
## 7 51 -5.830795
## 8 69 -5.475800
## 9 70 -4.144568
## 10 83 -3.328079
##
## $T3_WHS
## Row T3_WHS_mad
## 1 22 -3.083386
##
## $T2_CLS
## Row T2_CLS_mad
## 1 87 -3.480372
##
## $T3_CLS
## Row T3_CLS_mad
## 1 10 -3.06115
##
## $T2_PANAS_neg
## Row T2_PANAS_neg_mad
## 1 7 3.035208
## 2 44 5.058681
##
## $T2_memory.altruistic
## Row T2_memory.altruistic_mad
## 1 8 5.232062
## 2 18 3.080544
## 3 49 3.726311
## 4 56 8.054610
## 5 66 9.220174
## 6 77 3.873390
## 7 87 3.803484
## 8 90 4.521580
##
## $T2_memory.aggressive
## Row T2_memory.aggressive_mad
## 1 1 3.045000
## 2 2 5.060588
## 3 7 5.035664
## 4 30 4.291231
## 5 34 5.927863
## 6 48 3.351726
## 7 49 3.652856
## 8 68 4.764036
## 9 84 4.659505
data %>%
as.data.frame %>%
filter(T1_Group == "Reflection") %>%
find_mad(col.list, criteria = 3)## 33 outlier(s) based on 3 median absolute deviations for variable(s):
## T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T1_dehumanization, T2_dehumanization, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive
##
## The following participants were considered outliers for more than one variable:
##
## Row n
## 1 2 2
## 2 4 2
## 3 5 3
## 4 8 3
## 5 22 6
## 6 25 6
## 7 29 4
## 8 31 2
## 9 36 2
## 10 37 5
## 11 39 5
## 12 43 4
## 13 50 2
## 14 53 5
##
## Outliers per variable:
##
## $T1_blastduration
## Row T1_blastduration_mad
## 1 50 -3.01203
##
## $T1_blastintensity.duration
## Row T1_blastintensity.duration_mad
## 1 8 4.013772
## 2 41 3.607649
##
## $T2_blastintensity.duration
## Row T2_blastintensity.duration_mad
## 1 5 3.228501
## 2 8 7.975666
## 3 43 3.268278
##
## $T3_blastintensity.duration
## Row T3_blastintensity.duration_mad
## 1 5 3.777493
## 2 8 7.395048
## 3 36 3.236847
## 4 37 3.781831
## 5 42 5.320105
##
## $T2_attitude
## Row T2_attitude_mad
## 1 2 -3.489757
## 2 18 -11.524994
## 3 22 -5.190646
## 4 23 -3.357791
## 5 25 -5.410589
## 6 29 -6.407662
## 7 37 -5.615869
## 8 39 -3.123185
## 9 43 -3.357791
## 10 53 -4.105596
##
## $T3_attitude
## Row T3_attitude_mad
## 1 11 -7.503710
## 2 22 -6.070417
## 3 25 -5.817483
## 4 29 -6.306489
## 5 37 -5.261028
## 6 39 -5.328477
## 7 43 -4.974369
## 8 48 -3.878322
## 9 53 -7.874680
##
## $T1_dehumanization
## Row T1_dehumanization_mad
## 1 25 -10.791852
## 2 34 -3.736679
## 3 39 -4.910293
## 4 43 -5.935519
##
## $T2_dehumanization
## Row T2_dehumanization_mad
## 1 22 -4.903029
## 2 25 -4.903029
## 3 29 -4.928971
## 4 37 -4.474987
## 5 38 -4.124770
## 6 39 -4.215567
## 7 53 -3.891293
##
## $T3_dehumanization
## Row T3_dehumanization_mad
## 1 22 -3.839409
## 2 25 -3.320570
## 3 29 -3.663004
## 4 39 -3.901670
## 5 53 -3.548859
##
## $T1_WHS
## Row T1_WHS_mad
## 1 52 -3.597284
##
## $T2_SMS5
## Row T2_SMS5_mad
## 1 22 3.035208
## 2 50 3.709699
##
## $T2_PANAS_neg
## Row T2_PANAS_neg_mad
## 1 2 3.372454
## 2 12 3.372454
## 3 24 3.372454
## 4 53 4.046945
##
## $T2_Familiarity
## Row T2_Familiarity_mad
## 1 20 3.065867
## 2 36 3.065867
## 3 37 3.556406
##
## $T2_memory.altruistic
## Row T2_memory.altruistic_mad
## 1 3 3.003502
## 2 4 3.129885
## 3 30 3.723668
## 4 31 3.455823
## 5 49 3.000397
##
## $T2_memory.aggressive
## Row T2_memory.aggressive_mad
## 1 1 4.158773
## 2 4 6.946257
## 3 5 6.538862
## 4 15 3.383314
## 5 16 3.669782
## 6 22 11.217445
## 7 25 3.247418
## 8 31 5.680633
## 9 44 10.604004
data %>%
as.data.frame %>%
filter(T1_Group == "Meditation") %>%
find_mad(col.list, criteria = 3)## 36 outlier(s) based on 3 median absolute deviations for variable(s):
## T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T1_dehumanization, T2_dehumanization, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive
##
## The following participants were considered outliers for more than one variable:
##
## Row n
## 1 1 2
## 2 3 4
## 3 5 3
## 4 10 2
## 5 11 3
## 6 16 2
## 7 20 3
## 8 26 2
## 9 35 4
## 10 38 4
## 11 41 3
## 12 42 2
## 13 48 5
## 14 50 3
## 15 53 5
## 16 55 4
## 17 56 2
## 18 58 4
##
## Outliers per variable:
##
## $T1_blastintensity
## Row T1_blastintensity_mad
## 1 58 3.102657
##
## $T1_blastduration
## Row T1_blastduration_mad
## 1 19 -3.157726
## 2 27 -3.157726
## 3 38 -3.157726
## 4 58 3.268636
##
## $T1_blastintensity.duration
## Row T1_blastintensity.duration_mad
## 1 37 3.001374
## 2 53 3.041189
## 3 58 5.630293
##
## $T2_blastintensity.duration
## Row T2_blastintensity.duration_mad
## 1 6 3.004049
##
## $T1_attitude
## Row T1_attitude_mad
## 1 3 -6.936605
## 2 5 -5.715422
## 3 11 -4.437439
## 4 35 -5.218429
## 5 38 -3.869447
## 6 41 -5.289428
## 7 48 -5.630223
## 8 55 -4.153443
##
## $T2_attitude
## Row T2_attitude_mad
## 1 1 -3.019149
## 2 5 -6.552196
## 3 20 -6.552196
## 4 35 -6.231010
## 5 41 -5.733171
## 6 42 -9.619523
## 7 48 -6.423722
## 8 49 -3.982707
## 9 50 -4.769613
## 10 55 -5.990120
##
## $T3_attitude
## Row T3_attitude_mad
## 1 1 -3.495088
## 2 3 -10.460739
## 3 5 -10.362631
## 4 16 -6.610009
## 5 35 -10.362631
## 6 41 -5.825147
## 7 45 -3.568669
## 8 48 -9.234392
## 9 50 -4.083735
## 10 53 -6.094944
## 11 55 -10.411685
## 12 56 -3.372454
## 13 58 -10.779589
##
## $T1_dehumanization
## Row T1_dehumanization_mad
## 1 16 -3.912046
## 2 31 -8.129389
## 3 48 -6.027817
## 4 53 -6.084617
## 5 57 -3.528652
##
## $T2_dehumanization
## Row T2_dehumanization_mad
## 1 42 -9.011510
## 2 50 -3.443040
## 3 53 -3.097952
##
## $T3_dehumanization
## Row T3_dehumanization_mad
## 1 3 -5.039945
## 2 8 -5.658228
## 3 35 -7.381927
## 4 53 -5.377190
##
## $T1_CLS
## Row T1_CLS_mad
## 1 38 -3.460431
## 2 48 -4.398853
##
## $T3_CLS
## Row T3_CLS_mad
## 1 38 -3.059297
##
## $T2_SMS5
## Row T2_SMS5_mad
## 1 21 3.102657
##
## $T2_PANAS_neg
## Row T2_PANAS_neg_mad
## 1 13 3.709699
## 2 40 5.058681
## 3 54 3.372454
## 4 56 3.709699
##
## $T2_memory.altruistic
## Row T2_memory.altruistic_mad
## 1 3 3.917226
## 2 10 3.177224
## 3 11 7.148862
## 4 20 5.338032
## 5 23 3.700884
## 6 26 4.172745
## 7 28 3.457989
## 8 55 9.298788
##
## $T2_memory.aggressive
## Row T2_memory.aggressive_mad
## 1 2 7.108767
## 2 10 20.601909
## 3 11 3.306058
## 4 18 22.846607
## 5 20 14.980875
## 6 26 3.071248
## 7 34 3.698610
After our transformations, there are 12 outliers in the waitlist group, 3 in the reflection group, and 12 in the meditation group.
For multivariate outliers, it is recommended to use the Minimum Covariance Determinant, a robust version of the Mahalanobis distance (MCD, Leys et al., 2019).
Leys, C., Delacre, M., Mora, Y. L., Lakens, D., & Ley, C. (2019). How to classify, detect, and manage univariate and multivariate outliers, with emphasis on pre-registration. International Review of Social Psychology, 32(1).
data.na <- na.omit(data[col.list])
x <- check_outliers(data.na[-c(15:16)], method = "mcd", threshold = 500)
# We have to exclude two variables that are too large following the transformations
# Otherwise we get an error:
# Error in solve.default(cov, ...) :
# system is computationally singular: reciprocal condition number = 8.99059e-20
x## OK: No outliers detected.
## - Based on the following method and threshold: mcd (500).
## - For variables: T1_blastintensity, T1_blastduration, T1_blastintensity.duration, T2_blastintensity, T2_blastduration, T2_blastintensity.duration, T3_blastintensity, T3_blastduration, T3_blastintensity.duration, T1_BSCS, T1_BAQ, T1_attitude, T2_attitude, T3_attitude, T3_dehumanization, T1_NOBAGS, T1_WHS, T2_WHS, T3_WHS, T1_CLS, T2_CLS, T3_CLS, T2_SMS5, T2_PANAS_pos, T2_PANAS_neg, T2_Charity, T2_Familiarity, T2_memory.altruistic, T2_memory.aggressive
plot(x)There are 4 multivariate outliers according to the MCD method using an artificially high threshold. However, we did not mention in the preregistration that we would exclude multivariate outliers, so we will not at this time.
Visual assessment and the MAD method confirm we have some outlier values. We could ignore them but because they could have disproportionate influence on the models, one recommendation is to winsorize them by bringing the values at 3 SD. Instead of using the standard deviation around the mean, however, we use the absolute deviation around the median, as it is more robust to extreme observations. For a discussion, see:
Leys, C., Klein, O., Bernard, P., & Licata, L. (2013). Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median. Journal of Experimental Social Psychology, 49(4), 764–766. https://doi.org/10.1016/j.jesp.2013.03.013
Outliers are still present but were brought back within reasonable limits, where applicable.
We can now standardize our variables.
We are now ready to move to the analyses proper.
# If we decide to only center variables instead of standardizing them
# (as in the preregistration), then let's do this instead
data <- data %>%
mutate(across(all_of(gsub(".t.w.s", "", col.list)),
\(x) standardize(x, center = TRUE, scale = FALSE),
.names = "{col}"))In this section, we perform the planned contrasts.
# Specify the order of factor levels for "Group".
# Otherwise R will alphabetize them.
data$T1_Group <- factor(data$T1_Group, levels = c("Meditation", "Reflection", "Waitlist"))
# Define our dependent variables
DV <- data %>% select(T2_NOBAGS:T2_Charity) %>% names
# First column (which variable)
Variable <- rep(DV, each = 3)
# Second column (which comparison)
Comparison <- rep(c("MeditationvsCTR",
"ReflectionvsCTR",
"MeditationvsReflection"),
length(DV))
# 14 == number of DV
# Make list of all formulas
formulas <- c(
"T2_NOBAGS ~ T1_Group * T1_NOBAGS",
"T2_attitude ~ T1_Group * T1_attitude",
"T2_dehumanization ~ T1_Group * T1_dehumanization",
"T2_IAT ~ T1_Group * T1_IAT",
"T2_SMS5 ~ T1_Group",
"T2_PANAS_pos ~ T1_Group",
"T2_PANAS_neg ~ T1_Group",
"T2_blastintensity ~ T1_Group * T1_blastintensity",
"T2_blastduration ~ T1_Group * T1_blastduration",
"T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
"T2_memory.altruistic ~ T1_Group",
"T2_memory.aggressive ~ T1_Group",
"T2_WHS ~ T1_Group * T1_WHS",
"T2_CLS ~ T1_Group * T1_CLS",
"T2_Charity ~ T1_Group * T2_Familiarity"
)
# Make list of all models
models.list <- sapply(formulas, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)
## Attempt with nice_lm_contrasts
x <- nice_lm_contrasts(models.list, group = "T1_Group", data = data)
x.table <- nice_table(x, highlight = TRUE)
x.tableDependent Variable | Comparison | df | t | p | d | 95% CI |
|---|---|---|---|---|---|---|
T2_NOBAGS | Meditation - Waitlist | 197 | -1.47 | .144 | -0.30 | [-0.64, 0.03] |
Reflection - Waitlist | 197 | -1.40 | .162 | -0.21 | [-0.59, 0.13] | |
Meditation - Reflection | 197 | -0.03 | .979 | -0.09 | [-0.50, 0.29] | |
T2_attitude | Meditation - Waitlist | 197 | 0.82 | .411 | 0.17 | [-0.18, 0.51] |
Reflection - Waitlist | 197 | 1.25 | .213 | 0.08 | [-0.28, 0.43] | |
Meditation - Reflection | 197 | -0.41 | .684 | 0.09 | [-0.29, 0.50] | |
T2_dehumanization | Meditation - Waitlist | 195 | 0.82 | .415 | 0.13 | [-0.20, 0.39] |
Reflection - Waitlist | 195 | -0.65 | .518 | -0.15 | [-0.53, 0.22] | |
Meditation - Reflection | 195 | 1.31 | .191 | 0.28 | [-0.10, 0.65] | |
T2_IAT | Meditation - Waitlist | 197 | 0.51 | .611 | 0.08 | [-0.23, 0.39] |
Reflection - Waitlist | 197 | 1.14 | .256 | 0.38 | [0.01, 0.75] | |
Meditation - Reflection | 197 | -0.59 | .554 | -0.31 | [-0.68, 0.06] | |
T2_SMS5 | Meditation - Waitlist | 200 | -1.50 | .136 | -0.25 | [-0.55, 0.08] |
Reflection - Waitlist | 200 | -0.16 | .876 | -0.03 | [-0.38, 0.32] | |
Meditation - Reflection | 200 | -1.18 | .239 | -0.22 | [-0.59, 0.15] | |
T2_PANAS_pos | Meditation - Waitlist | 200 | 2.76 | .006** | 0.46 | [0.13, 0.82] |
Reflection - Waitlist | 200 | 1.17 | .244 | 0.20 | [-0.15, 0.53] | |
Meditation - Reflection | 200 | 1.38 | .170 | 0.26 | [-0.12, 0.65] | |
T2_PANAS_neg | Meditation - Waitlist | 200 | 1.60 | .111 | 0.27 | [-0.07, 0.62] |
Reflection - Waitlist | 200 | 0.44 | .662 | 0.08 | [-0.27, 0.40] | |
Meditation - Reflection | 200 | 1.01 | .312 | 0.19 | [-0.19, 0.57] | |
T2_blastintensity | Meditation - Waitlist | 197 | -1.35 | .179 | -0.03 | [-0.37, 0.29] |
Reflection - Waitlist | 197 | -4.16 | < .001*** | -0.19 | [-0.57, 0.19] | |
Meditation - Reflection | 197 | 2.62 | .009** | 0.16 | [-0.22, 0.57] | |
T2_blastduration | Meditation - Waitlist | 197 | -0.83 | .408 | -0.02 | [-0.34, 0.30] |
Reflection - Waitlist | 197 | -4.29 | < .001*** | -0.33 | [-0.68, 0.05] | |
Meditation - Reflection | 197 | 3.18 | .002** | 0.31 | [-0.08, 0.70] | |
T2_blastintensity.duration | Meditation - Waitlist | 197 | -1.06 | .291 | -0.04 | [-0.36, 0.27] |
Reflection - Waitlist | 197 | -3.20 | .002** | -0.14 | [-0.52, 0.25] | |
Meditation - Reflection | 197 | 1.99 | .048* | 0.10 | [-0.30, 0.51] | |
T2_memory.altruistic | Meditation - Waitlist | 200 | -0.52 | .603 | -0.09 | [-0.42, 0.30] |
Reflection - Waitlist | 200 | -2.44 | .015* | -0.42 | [-0.65, -0.15] | |
Meditation - Reflection | 200 | 1.76 | .080 | 0.33 | [0.06, 0.68] | |
T2_memory.aggressive | Meditation - Waitlist | 200 | 1.36 | .176 | 0.23 | [-0.16, 0.59] |
Reflection - Waitlist | 200 | 0.36 | .717 | 0.06 | [-0.21, 0.36] | |
Meditation - Reflection | 200 | 0.87 | .386 | 0.17 | [-0.27, 0.56] | |
T2_WHS | Meditation - Waitlist | 197 | 0.64 | .524 | 0.07 | [-0.25, 0.42] |
Reflection - Waitlist | 197 | 1.24 | .218 | -0.01 | [-0.35, 0.31] | |
Meditation - Reflection | 197 | -0.57 | .573 | 0.08 | [-0.31, 0.45] | |
T2_CLS | Meditation - Waitlist | 197 | 2.56 | .011* | 0.40 | [0.04, 0.72] |
Reflection - Waitlist | 197 | 3.46 | .001*** | 0.43 | [0.08, 0.77] | |
Meditation - Reflection | 197 | -0.88 | .382 | -0.03 | [-0.43, 0.37] | |
T2_Charity | Meditation - Waitlist | 189 | -1.55 | .124 | -0.19 | [-0.54, 0.15] |
Reflection - Waitlist | 189 | -0.93 | .355 | -0.09 | [-0.43, 0.25] | |
Meditation - Reflection | 189 | -0.52 | .601 | -0.10 | [-0.49, 0.34] |
# Make list of all formulas
formulas2 <- c(
"T3_NOBAGS ~ T1_Group * T1_NOBAGS",
"T3_attitude ~ T1_Group * T1_attitude",
"T3_dehumanization ~ T1_Group * T1_dehumanization",
"T3_IAT ~ T1_Group * T1_IAT",
"T3_blastintensity ~ T1_Group * T1_blastintensity",
"T3_blastduration ~ T1_Group * T1_blastduration",
"T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
"T3_WHS ~ T1_Group * T1_WHS",
"T3_CLS ~ T1_Group * T1_CLS"
)
# Make list of all models
models.list2 <- sapply(formulas2, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)
## Attempt with nice_lm_contrasts
x2 <- nice_lm_contrasts(models.list2, group = "T1_Group", data = data)
x.table2 <- nice_table(x2, highlight = TRUE)
x.table2Dependent Variable | Comparison | df | t | p | d | 95% CI |
|---|---|---|---|---|---|---|
T3_NOBAGS | Meditation - Waitlist | 197 | -1.05 | .297 | -0.24 | [-0.57, 0.06] |
Reflection - Waitlist | 197 | 0.34 | .735 | 0.05 | [-0.31, 0.41] | |
Meditation - Reflection | 197 | -1.23 | .220 | -0.29 | [-0.70, 0.09] | |
T3_attitude | Meditation - Waitlist | 197 | 0.88 | .379 | 0.19 | [-0.16, 0.51] |
Reflection - Waitlist | 197 | 1.44 | .150 | 0.07 | [-0.26, 0.41] | |
Meditation - Reflection | 197 | -0.53 | .594 | 0.11 | [-0.28, 0.48] | |
T3_dehumanization | Meditation - Waitlist | 196 | 1.23 | .221 | 0.16 | [-0.16, 0.42] |
Reflection - Waitlist | 196 | -0.82 | .415 | -0.17 | [-0.55, 0.22] | |
Meditation - Reflection | 196 | 1.83 | .069 | 0.33 | [-0.05, 0.72] | |
T3_IAT | Meditation - Waitlist | 197 | 1.30 | .196 | 0.20 | [-0.09, 0.55] |
Reflection - Waitlist | 197 | 1.15 | .250 | 0.36 | [0.02, 0.72] | |
Meditation - Reflection | 197 | 0.08 | .938 | -0.16 | [-0.55, 0.22] | |
T3_blastintensity | Meditation - Waitlist | 197 | -1.79 | .076 | -0.11 | [-0.44, 0.22] |
Reflection - Waitlist | 197 | -2.82 | .005** | -0.11 | [-0.46, 0.25] | |
Meditation - Reflection | 197 | 1.02 | .310 | -0.01 | [-0.38, 0.40] | |
T3_blastduration | Meditation - Waitlist | 197 | -0.94 | .348 | -0.06 | [-0.36, 0.27] |
Reflection - Waitlist | 197 | -2.69 | .008** | -0.21 | [-0.58, 0.13] | |
Meditation - Reflection | 197 | 1.62 | .106 | 0.15 | [-0.25, 0.51] | |
T3_blastintensity.duration | Meditation - Waitlist | 197 | -1.16 | .247 | -0.08 | [-0.40, 0.23] |
Reflection - Waitlist | 197 | -2.13 | .034* | -0.06 | [-0.41, 0.32] | |
Meditation - Reflection | 197 | 0.93 | .354 | -0.02 | [-0.42, 0.35] | |
T3_WHS | Meditation - Waitlist | 197 | 1.68 | .095 | 0.17 | [-0.15, 0.49] |
Reflection - Waitlist | 197 | 2.74 | .007** | 0.11 | [-0.23, 0.45] | |
Meditation - Reflection | 197 | -1.03 | .307 | 0.06 | [-0.30, 0.44] | |
T3_CLS | Meditation - Waitlist | 196 | 1.99 | .048* | 0.33 | [-0.01, 0.66] |
Reflection - Waitlist | 196 | 3.03 | .003** | 0.37 | [-0.02, 0.69] | |
Meditation - Reflection | 196 | -0.99 | .322 | -0.04 | [-0.43, 0.33] |
In this section, we generate various plots: violin plots, marginal means plots, and lighthouse plots, to better illustrate the pairwise contrasts (only for significant contrasts from the previous step).
nice_violin(data,
group = "T1_Group",
response = "T2_attitude",
obs = TRUE,
signif_annotation = c("*"),
signif_yposition = 1.5,
signif_xmin = 2,
signif_xmax = 3)The first plot below represents the adjusted (marginal) means, after taking into account covariates such as Time 1 measurements.
means <- estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Attitude", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(models.list$`T2_attitude ~ T1_Group * T1_attitude`)
plot(contrasts, estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)) +
ylab(paste("Adjusted", "Attitude", "Mean")) +
theme_modern()The plot above is a lighthouse plot. These plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T2_PANAS_pos",
obs = TRUE,
comp1 = 1,
comp2 = 3,
has.d = TRUE,
d.x = 1.75,
d.y = 2)means <- estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Positive Affect", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(models.list$`T2_PANAS_pos ~ T1_Group`)
plot(contrasts, estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)) +
ylab(paste("Adjusted", "Positive Affect", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T2_blastintensity",
obs = TRUE,
signif_annotation = c("**", "***"),
signif_yposition = c(3.5, 4),
signif_xmin = c(1, 2),
signif_xmax = c(2, 3))means <- estimate_means(models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T2_blastduration",
obs = TRUE,
signif_annotation = c("**", "***"),
signif_yposition = c(3.5, 4),
signif_xmin = c(1, 2),
signif_xmax = c(2, 3))means <- estimate_means(models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T2_blastintensity.duration",
obs = TRUE,
signif_annotation = c("*", "**"),
signif_yposition = c(3.5, 4),
signif_xmin = c(1, 2),
signif_xmax = c(2, 3))means <- estimate_means(
models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T2_memory.altruistic",
obs = TRUE,
signif_annotation = c("*"),
signif_yposition = 3.5,
signif_xmin = 2,
signif_xmax = 3)means <- estimate_means(
models.list$`T2_memory.altruistic ~ T1_Group`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list$`T2_memory.altruistic ~ T1_Group`)
plot(contrasts, estimate_means(
models.list$`T2_memory.altruistic ~ T1_Group`)) +
ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T2_CLS",
obs = TRUE,
signif_annotation = c("*", "***"),
signif_yposition = c(3.2, 2.5),
signif_xmin = c(1, 2),
signif_xmax = c(3, 3))means <- estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Compassion", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(models.list$`T2_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T3_blastintensity",
obs = TRUE,
signif_annotation = c("**"),
signif_yposition = 2.5,
signif_xmin = 2,
signif_xmax = 3)means <- estimate_means(models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T3_blastduration",
obs = TRUE,
signif_annotation = c("**"),
signif_yposition = 3,
signif_xmin = 2,
signif_xmax = 3)means <- estimate_means(models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()nice_violin(data,
group = "T1_Group",
response = "T3_blastintensity.duration",
obs = TRUE,
signif_annotation = c("*"),
signif_yposition = 3,
signif_xmin = 2,
signif_xmax = 3)means <- estimate_means(
models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T3_WHS",
obs = TRUE,
signif_annotation = c("**"),
signif_yposition = 2.5,
signif_xmin = 2,
signif_xmax = 3)means <- estimate_means(models.list2$`T3_WHS ~ T1_Group * T1_WHS`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(
models.list2$`T3_WHS ~ T1_Group * T1_WHS`)
plot(contrasts, estimate_means(
models.list2$`T3_WHS ~ T1_Group * T1_WHS`)) +
ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
nice_violin(data,
group = "T1_Group",
response = "T3_CLS",
obs = TRUE,
signif_annotation = c("*", "**"),
signif_yposition = c(3.2, 2.7),
signif_xmin = c(1, 2),
signif_xmax = c(3, 3))means <- estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Compassion", "Mean")) +
theme_modern()contrasts <- estimate_contrasts(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()Lighthouse plots represent the estimated means and their CI range (in black), while the grey areas show the CI range of the difference (as compared to the point estimate).
When conducting a Randomized-Controlled Trial (RCT), some participants are randomly assigned to a treatment condition, and others to a control group. However, not everyone assigned to the treatment might follow the treatment protocol (called “treatment compliance”).
According to Sagarin et al. (2014), one sensible approach to address this problem is using the complier average causal effect (CACE), also sometimes known as Local average treatment effect (LATE). According to Wikipedia, it is “the treatment effect for the subset of the sample that takes the treatment if and only if they were assigned to the treatment, otherwise known as the compliers.” In other words, it will be useful if a proportion of your participants assigned to the treatment group did not follow the treatment protocol.
In the following figures, the P > in the compliance
column represents the compliance Percentage. The numbers on the right
represent Hedge’s g (analogous to Cohen’s d) and its
95% confidence interval.
data3 <- data2 %>%
mutate(part.percent = ifelse(T1_Group == "Waitlist", 0, part.percent),
part.percent = part.percent * 100,
T2_memory.altruistic = T2_memory.altruistic * -1,
T1_GroupReflection = as.numeric(T1_Group == "Reflection"),
T1_GroupMeditation = as.numeric(T1_Group == "Meditation"),
across(contains("blast"), \(x) x * -1)) %>%
as.data.frame()
caceOutput <- caceSRTBoot(
T2_attitude ~ T1_attitude + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T2_PANAS_pos ~ T1_GroupMeditation,
intervention = "T1_GroupMeditation",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Positive Affect (Meditation VS Waitlist)")caceOutput <- caceSRTBoot(
T2_blastintensity ~ T1_blastintensity + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T2_blastduration ~ T1_blastduration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T2_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T2_memory.altruistic ~ T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "ms to remember altruistic event (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T2_CLS ~ T1_CLS + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T2_CLS ~ T1_CLS + T1_GroupMeditation,
intervention = "T1_GroupMeditation",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")caceOutput <- caceSRTBoot(
T3_attitude ~ T1_attitude + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T3_blastintensity ~ T1_blastintensity + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T3_blastduration ~ T1_blastduration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T3_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T3_WHS ~ T1_WHS + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Willingness to Help (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T3_CLS ~ T1_CLS + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")caceOutput <- caceSRTBoot(
T3_CLS ~ T1_CLS + T1_GroupMeditation,
intervention = "T1_GroupMeditation",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")What we see is general trend toward larger effects with higher compliance, although this trend does not appear exceptionally strong.
It seems like only IAT (implicit aggression) and not NOBAGS (explicit attitude toward aggression) is a significant moderator.
# CREATE OUR DUMMY VARIABLES FOR T1_Group!
data$T1_GroupReflection <- as.numeric(data$T1_Group == "Reflection")
data$T1_GroupMeditation <- as.numeric(data$T1_Group == "Meditation")
################################################
################################################
# T2_blastintensity.duration - IAT
T2_blastintensity.duration.IAT <- lm(
T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition +
T1_GroupReflection:T2_IAT:T2_Condition + T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_blastintensity.duration.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_blastintensity.duration.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_blastintensity.duration | T1_GroupReflection | 191 | -405.54 | -0.31 | .760 | .00 | [0.00, 0.01] |
T1_GroupMeditation | 191 | -54.71 | -0.03 | .974 | .00 | [0.00, 0.00] | |
T2_IAT | 191 | -155.01 | -0.11 | .915 | .00 | [0.00, 0.00] | |
T2_ConditionDepleted | 191 | 2,409.51 | 1.68 | .096 | .01 | [0.00, 0.04] | |
T1_GroupReflection × T2_IAT | 191 | 215.91 | 0.10 | .917 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT | 191 | -69.14 | -0.03 | .979 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | -2,181.85 | -0.85 | .398 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -5,380.99 | -2.14 | .033* | .02 | [0.00, 0.06] | |
T2_IAT × T2_ConditionDepleted | 191 | 2,746.44 | 1.41 | .161 | .01 | [0.00, 0.04] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | -3,823.21 | -1.06 | .290 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | -7,357.76 | -2.06 | .041* | .02 | [0.00, 0.06] |
# Make interaction plot
interact_plot(T2_blastintensity.duration.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Blast Intensity * Duration (Taylor)",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))What we see here is that for the waitlist group, there is no interaction between being depleted and implicit aggression, whereas for the meditation group, people become less aggressive after being depleted, the higher their implicit aggression.
# T2_blastintensity - IAT
T2_blastintensity.IAT <- lm(
T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_blastintensity.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_blastintensity.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_blastintensity | T1_GroupReflection | 191 | -0.59 | -0.67 | .502 | .00 | [0.00, 0.02] |
T1_GroupMeditation | 191 | -0.19 | -0.17 | .864 | .00 | [0.00, 0.00] | |
T2_IAT | 191 | -0.09 | -0.09 | .928 | .00 | [0.00, 0.00] | |
T2_ConditionDepleted | 191 | 0.84 | 0.88 | .381 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_IAT | 191 | -0.07 | -0.05 | .961 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT | 191 | -0.07 | -0.04 | .966 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | -0.79 | -0.46 | .646 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -2.29 | -1.38 | .171 | .01 | [0.00, 0.04] | |
T2_IAT × T2_ConditionDepleted | 191 | 1.08 | 0.83 | .406 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | -1.88 | -0.79 | .432 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | -3.49 | -1.47 | .142 | .01 | [0.00, 0.04] |
# T2_blastduration - IAT
T2_blastduration.IAT <- lm(
T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_blastduration.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_blastduration.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_blastduration | T1_GroupReflection | 191 | -151.94 | -0.84 | .403 | .00 | [0.00, 0.02] |
T1_GroupMeditation | 191 | 133.46 | 0.58 | .560 | .00 | [0.00, 0.01] | |
T2_IAT | 191 | -92.44 | -0.47 | .641 | .00 | [0.00, 0.01] | |
T2_ConditionDepleted | 191 | 271.17 | 1.38 | .169 | .01 | [0.00, 0.04] | |
T1_GroupReflection × T2_IAT | 191 | 8.01 | 0.03 | .977 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT | 191 | 196.17 | 0.55 | .582 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | -231.85 | -0.66 | .511 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -706.58 | -2.06 | .041* | .02 | [0.00, 0.06] | |
T2_IAT × T2_ConditionDepleted | 191 | 326.41 | 1.22 | .223 | .01 | [0.00, 0.03] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | -471.58 | -0.96 | .339 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | -962.19 | -1.97 | .050 | .02 | [0.00, 0.06] |
# Make interaction plot
interact_plot(T2_blastduration.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Blast Duration (Taylor)",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))What we see here is that for the waitlist group, there is no interaction between being depleted and implicit aggression, whereas for the meditation group, people become less aggressive after being depleted, the higher their implicit aggression. However, that variable (blastduration alone) was not in the preregistration, so we might not report this finding.
# T2_Charity - IAT
T2_Charity.IAT <- lm(
T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition + T2_Familiarity,
data = data)
check_model(T2_Charity.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_Charity.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_Charity | T1_GroupReflection | 182 | 0.91 | 0.71 | .480 | .00 | [0.00, 0.02] |
T1_GroupMeditation | 182 | -0.72 | -0.43 | .670 | .00 | [0.00, 0.01] | |
T2_IAT | 182 | -0.47 | -0.34 | .732 | .00 | [0.00, 0.01] | |
T2_ConditionDepleted | 182 | 1.04 | 0.76 | .451 | .00 | [0.00, 0.02] | |
T2_Familiarity | 182 | 1.04 | 3.56 | < .001*** | .06 | [0.00, 0.13] | |
T1_GroupReflection × T2_IAT | 182 | 1.55 | 0.78 | .435 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_IAT | 182 | -0.26 | -0.10 | .921 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_ConditionDepleted | 182 | -2.03 | -0.79 | .430 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_ConditionDepleted | 182 | -0.31 | -0.13 | .900 | .00 | [0.00, 0.00] | |
T2_IAT × T2_ConditionDepleted | 182 | 1.32 | 0.70 | .487 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 182 | -0.78 | -0.22 | .825 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 182 | 0.12 | 0.03 | .974 | .00 | [0.00, 0.00] |
# Compassion - IAT
T2_CLS.IAT <- lm(
T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_CLS.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_CLS.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_CLS | T1_GroupReflection | 191 | 1.43 | 3.77 | < .001*** | .06 | [0.00, 0.12] |
T1_GroupMeditation | 191 | 1.42 | 2.97 | .003** | .04 | [0.00, 0.09] | |
T2_IAT | 191 | -1.47 | -3.54 | < .001*** | .06 | [0.00, 0.11] | |
T2_ConditionDepleted | 191 | 1.15 | 2.80 | .006** | .03 | [0.00, 0.08] | |
T1_GroupReflection × T2_IAT | 191 | 1.69 | 2.86 | .005** | .04 | [0.00, 0.08] | |
T1_GroupMeditation × T2_IAT | 191 | 1.77 | 2.37 | .019* | .02 | [0.00, 0.06] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 0.09 | 0.12 | .907 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -1.87 | -2.61 | .010** | .03 | [0.00, 0.07] | |
T2_IAT × T2_ConditionDepleted | 191 | 1.96 | 3.51 | .001*** | .05 | [0.00, 0.11] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | -0.03 | -0.03 | .979 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | -3.07 | -3.00 | .003** | .04 | [0.00, 0.09] |
# Make interaction plot
interact_plot(T2_CLS.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Compassionate Love",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))What we see here is that for the waitlist group, the effect of implicit aggression clearly depends on depletion: implicit aggression leads to lower compassion in the control condition (expected), but to higher compassion in the depletion group (unexpected). Whereas for the meditation group, it seems the depletion has little effect on compassion (but the trend is reversed relative to the waitlist group).
# T2_WHS - IAT
T2_WHS.IAT <- lm(
T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + T2_Condition +
T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_WHS.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_WHS.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_WHS | T1_GroupReflection | 191 | 0.53 | 1.43 | .154 | .01 | [0.00, 0.04] |
T1_GroupMeditation | 191 | 0.69 | 1.48 | .141 | .01 | [0.00, 0.04] | |
T2_IAT | 191 | -0.73 | -1.81 | .072 | .02 | [0.00, 0.05] | |
T2_ConditionDepleted | 191 | 0.50 | 1.25 | .213 | .01 | [0.00, 0.03] | |
T1_GroupReflection × T2_IAT | 191 | 0.90 | 1.57 | .119 | .01 | [0.00, 0.04] | |
T1_GroupMeditation × T2_IAT | 191 | 0.95 | 1.31 | .193 | .01 | [0.00, 0.03] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 0.38 | 0.53 | .595 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -0.88 | -1.26 | .209 | .01 | [0.00, 0.03] | |
T2_IAT × T2_ConditionDepleted | 191 | 0.85 | 1.55 | .123 | .01 | [0.00, 0.04] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | 0.49 | 0.49 | .624 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | -1.25 | -1.25 | .212 | .01 | [0.00, 0.03] |
# T2_memory.altruistic - IAT
T2_memory.altruistic.IAT <- lm(
T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_memory.altruistic.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_memory.altruistic.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_memory.altruistic | T1_GroupReflection | 191 | -3,870.57 | -2.08 | .039* | .02 | [0.00, 0.06] |
T1_GroupMeditation | 191 | -4,129.04 | -1.76 | .081 | .02 | [0.00, 0.05] | |
T2_IAT | 191 | 2,529.17 | 1.24 | .216 | .01 | [0.00, 0.03] | |
T2_ConditionDepleted | 191 | -4,289.34 | -2.12 | .035* | .02 | [0.00, 0.06] | |
T1_GroupReflection × T2_IAT | 191 | -2,907.13 | -1.00 | .318 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_IAT | 191 | -5,638.28 | -1.54 | .126 | .01 | [0.00, 0.04] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 3,963.65 | 1.10 | .275 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | 4,949.85 | 1.40 | .162 | .01 | [0.00, 0.04] | |
T2_IAT × T2_ConditionDepleted | 191 | -5,349.23 | -1.95 | .053 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | 5,764.69 | 1.14 | .256 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | 6,772.42 | 1.35 | .179 | .01 | [0.00, 0.03] |
# Make interaction plot
interact_plot(T2_memory.altruistic.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Reaction Time (Altruistic Memory)",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))What we see here is that for the waitlist group, the depletion seems to have little effect on reaction time to remember an altruistic event. Whereas for the meditation group, higher implicit aggression relates to shorter reaction time (unexpected), unless they are depleted (but the trend is reversed relative to the waitlist group). This result is somewhat unexpected, but it could be that in their case, the heart takes over the mind.
# T2_memory.aggressive - IAT
T2_memory.aggressive.IAT <- lm(
T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_memory.aggressive.IAT)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_memory.aggressive.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_memory.aggressive | T1_GroupReflection | 191 | 2,692.41 | 0.61 | .545 | .00 | [0.00, 0.01] |
T1_GroupMeditation | 191 | 4,552.33 | 0.81 | .417 | .00 | [0.00, 0.02] | |
T2_IAT | 191 | -2,085.63 | -0.43 | .668 | .00 | [0.00, 0.01] | |
T2_ConditionDepleted | 191 | 1,809.33 | 0.38 | .707 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_IAT | 191 | 1,263.35 | 0.18 | .855 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT | 191 | 5,654.74 | 0.65 | .518 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | -2,341.73 | -0.27 | .786 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -7,225.39 | -0.86 | .391 | .00 | [0.00, 0.02] | |
T2_IAT × T2_ConditionDepleted | 191 | -87.36 | -0.01 | .989 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_IAT × T2_ConditionDepleted | 191 | 224.90 | 0.02 | .985 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_IAT × T2_ConditionDepleted | 191 | -14,323.47 | -1.20 | .233 | .01 | [0.00, 0.03] |
# T2_blastintensity.duration - NOBAGS
T2_blastintensity.duration.NOBAGS <- lm(
T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_blastintensity.duration.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_blastintensity.duration.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_blastintensity.duration | T1_GroupReflection | 191 | -3,581.79 | -1.16 | .246 | .01 | [0.00, 0.03] |
T1_GroupMeditation | 191 | -974.72 | -0.29 | .770 | .00 | [0.00, 0.01] | |
T2_NOBAGS | 191 | -733.10 | -0.70 | .488 | .00 | [0.00, 0.02] | |
T2_ConditionDepleted | 191 | -5,165.29 | -1.60 | .112 | .01 | [0.00, 0.04] | |
T1_GroupReflection × T2_NOBAGS | 191 | 1,525.17 | 1.02 | .309 | .01 | [0.00, 0.02] | |
T1_GroupMeditation × T2_NOBAGS | 191 | 451.50 | 0.28 | .781 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 8,010.38 | 1.63 | .104 | .01 | [0.00, 0.04] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | 1,410.90 | 0.29 | .769 | .00 | [0.00, 0.01] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | 2,862.52 | 1.82 | .070 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | -3,841.98 | -1.56 | .121 | .01 | [0.00, 0.04] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | -613.81 | -0.25 | .800 | .00 | [0.00, 0.01] |
# T2_blastintensity - NOBAGS
T2_blastintensity.NOBAGS <- lm(
T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS +
T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_blastintensity.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_blastintensity.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_blastintensity | T1_GroupReflection | 191 | -2.72 | -1.35 | .177 | .01 | [0.00, 0.03] |
T1_GroupMeditation | 191 | -0.47 | -0.21 | .831 | .00 | [0.00, 0.00] | |
T2_NOBAGS | 191 | -0.55 | -0.80 | .422 | .00 | [0.00, 0.02] | |
T2_ConditionDepleted | 191 | -3.82 | -1.81 | .072 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_NOBAGS | 191 | 1.07 | 1.10 | .274 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_NOBAGS | 191 | 0.14 | 0.13 | .899 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 4.54 | 1.42 | .158 | .01 | [0.00, 0.04] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | 0.70 | 0.22 | .824 | .00 | [0.00, 0.00] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | 1.96 | 1.91 | .057 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | -1.97 | -1.22 | .223 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | -0.11 | -0.07 | .947 | .00 | [0.00, 0.00] |
# T2_blastduration - NOBAGS
T2_blastduration.NOBAGS <- lm(
T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS +
T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_blastduration.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_blastduration.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_blastduration | T1_GroupReflection | 191 | -491.94 | -1.20 | .233 | .01 | [0.00, 0.03] |
T1_GroupMeditation | 191 | 26.08 | 0.06 | .953 | .00 | [0.00, 0.00] | |
T2_NOBAGS | 191 | -119.38 | -0.85 | .397 | .00 | [0.00, 0.02] | |
T2_ConditionDepleted | 191 | -822.26 | -1.90 | .058 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_NOBAGS | 191 | 157.43 | 0.79 | .431 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_NOBAGS | 191 | -8.25 | -0.04 | .970 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 738.07 | 1.13 | .261 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -243.02 | -0.38 | .705 | .00 | [0.00, 0.01] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | 438.08 | 2.09 | .038* | .02 | [0.00, 0.06] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | -302.62 | -0.92 | .359 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | 143.03 | 0.44 | .659 | .00 | [0.00, 0.01] |
# Make interaction plot
interact_plot(T2_blastduration.NOBAGS, pred = "T2_NOBAGS", modx = "T2_Condition",
interval = TRUE,
x.label = "Normative beliefs about aggression (NOBAGS)",
y.label = "Blast Duration (Taylor)",
legend.main = "Condition") +
theme(text = element_text(size = 25))What we see here is that there seems to be no relation between attitude toward aggression on behavioural aggression in the control condition, whereas when depleted, attitude toward aggression relates to higher behavioural aggression. Although this result is theoretically consistent with the literature, it is likely a false positive given our high number of tests, the fact that this is the only variable that NOBAGS moderates, and that the p value is relatively close to 0.5. Furthermore that variable (blastduration alone) was not in the preregistration, so we might not report this finding.
# T2_Charity - NOBAGS
T2_Charity.NOBAGS <- lm(
T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition + T2_Familiarity,
data = data)
check_model(T2_Charity.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_Charity.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_Charity | T1_GroupReflection | 182 | 1.16 | 0.39 | .695 | .00 | [0.00, 0.01] |
T1_GroupMeditation | 182 | 1.13 | 0.36 | .721 | .00 | [0.00, 0.01] | |
T2_NOBAGS | 182 | 0.62 | 0.62 | .534 | .00 | [0.00, 0.01] | |
T2_ConditionDepleted | 182 | 2.98 | 0.96 | .339 | .00 | [0.00, 0.02] | |
T2_Familiarity | 182 | 1.17 | 3.93 | < .001*** | .08 | [0.01, 0.15] | |
T1_GroupReflection × T2_NOBAGS | 182 | -0.51 | -0.35 | .725 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_NOBAGS | 182 | -0.86 | -0.56 | .577 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_ConditionDepleted | 182 | -8.02 | -1.69 | .092 | .01 | [0.00, 0.05] | |
T1_GroupMeditation × T2_ConditionDepleted | 182 | -5.03 | -1.09 | .275 | .01 | [0.00, 0.03] | |
T2_NOBAGS × T2_ConditionDepleted | 182 | -1.38 | -0.92 | .359 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 182 | 3.29 | 1.37 | .174 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 182 | 2.40 | 1.03 | .303 | .01 | [0.00, 0.02] |
# Compassion - NOBAGS
T2_CLS.NOBAGS <- lm(
T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_CLS.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_CLS.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_CLS | T1_GroupReflection | 191 | 0.13 | 0.14 | .889 | .00 | [0.00, 0.00] |
T1_GroupMeditation | 191 | 1.51 | 1.54 | .125 | .01 | [0.00, 0.04] | |
T2_NOBAGS | 191 | -0.24 | -0.76 | .448 | .00 | [0.00, 0.02] | |
T2_ConditionDepleted | 191 | 0.63 | 0.67 | .505 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_NOBAGS | 191 | 0.14 | 0.31 | .757 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_NOBAGS | 191 | -0.61 | -1.27 | .204 | .01 | [0.00, 0.03] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | -0.27 | -0.19 | .850 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -1.68 | -1.19 | .235 | .01 | [0.00, 0.03] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | -0.38 | -0.82 | .414 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | 0.10 | 0.14 | .892 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | 0.88 | 1.24 | .215 | .01 | [0.00, 0.03] |
# T2_WHS - NOBAGS
T2_WHS.NOBAGS <- lm(
T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + T2_Condition +
T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_WHS.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_WHS.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_WHS | T1_GroupReflection | 191 | -1.04 | -1.22 | .222 | .01 | [0.00, 0.03] |
T1_GroupMeditation | 191 | 0.92 | 1.01 | .316 | .00 | [0.00, 0.02] | |
T2_NOBAGS | 191 | -0.47 | -1.63 | .104 | .01 | [0.00, 0.04] | |
T2_ConditionDepleted | 191 | -1.67 | -1.87 | .063 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_NOBAGS | 191 | 0.51 | 1.23 | .219 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_NOBAGS | 191 | -0.43 | -0.97 | .335 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 1.90 | 1.41 | .161 | .01 | [0.00, 0.04] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -0.54 | -0.41 | .683 | .00 | [0.00, 0.01] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | 0.81 | 1.87 | .063 | .02 | [0.00, 0.05] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | -0.97 | -1.43 | .154 | .01 | [0.00, 0.04] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | 0.26 | 0.40 | .693 | .00 | [0.00, 0.01] |
# T2_memory.altruistic - NOBAGS
T2_memory.altruistic.NOBAGS <- lm(
T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS +
T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_memory.altruistic.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_memory.altruistic.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_memory.altruistic | T1_GroupReflection | 191 | -1,307.87 | -0.30 | .761 | .00 | [0.00, 0.01] |
T1_GroupMeditation | 191 | 1,767.77 | 0.38 | .705 | .00 | [0.00, 0.01] | |
T2_NOBAGS | 191 | 1,260.19 | 0.86 | .393 | .00 | [0.00, 0.02] | |
T2_ConditionDepleted | 191 | -3,774.82 | -0.84 | .404 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_NOBAGS | 191 | -360.53 | -0.17 | .863 | .00 | [0.00, 0.00] | |
T1_GroupMeditation × T2_NOBAGS | 191 | -1,283.78 | -0.57 | .572 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | 7,532.92 | 1.10 | .273 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | 1,944.98 | 0.29 | .772 | .00 | [0.00, 0.01] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | 1,527.83 | 0.70 | .487 | .00 | [0.00, 0.02] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | -3,795.29 | -1.10 | .272 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | -415.11 | -0.12 | .903 | .00 | [0.00, 0.00] |
# T2_memory.aggressive - NOBAGS
T2_memory.aggressive.NOBAGS <- lm(
T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_memory.aggressive.NOBAGS)## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
T2_memory.aggressive.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)Dependent Variable | Predictor | df | b | t | p | sr2 | 95% CI |
|---|---|---|---|---|---|---|---|
T2_memory.aggressive | T1_GroupReflection | 191 | 8,077.14 | 0.79 | .429 | .00 | [0.00, 0.02] |
T1_GroupMeditation | 191 | -5,022.93 | -0.46 | .650 | .00 | [0.00, 0.01] | |
T2_NOBAGS | 191 | -1,512.22 | -0.43 | .665 | .00 | [0.00, 0.01] | |
T2_ConditionDepleted | 191 | 1,092.80 | 0.10 | .919 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_NOBAGS | 191 | -3,290.98 | -0.67 | .506 | .00 | [0.00, 0.01] | |
T1_GroupMeditation × T2_NOBAGS | 191 | 3,160.31 | 0.59 | .557 | .00 | [0.00, 0.01] | |
T1_GroupReflection × T2_ConditionDepleted | 191 | -17,673.12 | -1.09 | .278 | .01 | [0.00, 0.03] | |
T1_GroupMeditation × T2_ConditionDepleted | 191 | -4,542.27 | -0.29 | .775 | .00 | [0.00, 0.01] | |
T2_NOBAGS × T2_ConditionDepleted | 191 | 395.27 | 0.08 | .939 | .00 | [0.00, 0.00] | |
T1_GroupReflection × T2_NOBAGS × T2_ConditionDepleted | 191 | 8,078.65 | 0.99 | .323 | .00 | [0.00, 0.02] | |
T1_GroupMeditation × T2_NOBAGS × T2_ConditionDepleted | 191 | 4,027.85 | 0.50 | .616 | .00 | [0.00, 0.01] |
data %>%
select(T1_attitude_1:T1_attitude_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Positive Explicit Attitude",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")data %>%
select(T1_dehumanization_1:T1_dehumanization_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Humanization",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")data %>%
select(T2_attitude_1:T2_attitude_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Positive Explicit Attitude",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")data %>%
select(T2_dehumanization_1:T2_dehumanization_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Humanization",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")data %>%
select(contains("charity") & ends_with("1_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "donation") %>%
mutate(charity = factor(charity, labels = charities)) %>%
nice_violin(group = "charity",
response = "donation",
ytitle = "Amount Donated",
CIcap.width = 0.5,
obs = "jitter",
order.groups = "increasing",
xlabels.angle = 75)data %>%
select(contains("charity") & ends_with("2_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "familiarity") %>%
mutate(charity = factor(charity, labels = charities)) %>%
nice_violin(group = "charity",
response = "familiarity",
ytitle = "Familiarity with Charity",
CIcap.width = 0.5,
obs = "jitter",
order.groups = "increasing",
xlabels.angle = 75)data %>%
nice_scatter(predictor = "T2_Familiarity",
response = "T2_Charity",
ytitle = "Donation Amount",
xtitle = "Familiarity with Charity",
has.jitter = TRUE,
has.legend = TRUE,
has.r = TRUE,
has.p = TRUE)data %>%
select(contains("charity") & ends_with("1_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "donation") %>%
mutate(charity = factor(charity, labels = charities),
region = case_match(
charity,
charities[1:6] ~ regions[1],
charities[7:12] ~ regions[2],
charities[13:18] ~ regions[3],
charities[19:24] ~ regions[4]
)) %>%
nice_violin(group = "region",
response = "donation",
ytitle = "Amount Donated",
obs = "jitter",
order.groups = "increasing")data %>%
select(contains("charity") & ends_with("2_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "familiarity") %>%
mutate(charity = factor(charity, labels = charities),
region = case_match(
charity,
charities[1:6] ~ regions[1],
charities[7:12] ~ regions[2],
charities[13:18] ~ regions[3],
charities[19:24] ~ regions[4]
)) %>%
nice_violin(group = "region",
response = "familiarity",
ytitle = "Familiarity with Charity",
obs = "jitter",
order.groups = "increasing")data %>%
select(T3_attitude_1:T3_attitude_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Positive Explicit Attitude",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")data %>%
select(T3_dehumanization_1:T3_dehumanization_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Humanization",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")In this report, we aimed to compare two types of loving-kindness meditation, one more embodied, based on meditation, and one more cognitive, based on intellectual reflection, to a waitlist control group. We compared those groups on several variables relating to prosociality (i.e., on affect, attitude and behaviour). Groups were measured three times: Week 0 (T1), Week 6 (T2), and Week 13 (T3) so it was able to compare for baseline but also see how robust the effects, if any, are through times. We were also interested in assessing whether these effects depend on other personality variables (i.e., moderators).
Our contrasts analyses first revealed group differences at Time 2. Both the meditation and reflection groups showed moderately more compassionate love than the waitlist group, but only the meditation group showed moderately more positive affect than the waitlist group. However, the reflection group showed a little more positive explicit attitudes toward various social groups, as well as moderately shorter reaction times to remember an altruistic event than the waitlist group (suggesting that altruism was more cognitively accessible to them). Furthermore, the reflection group showed a little lower behavioural aggression (blast intensity, blast duration, and blast intensity * duration) than both the waitlist group and the meditation group.
Our contrasts analyses also revealed group differences at Time 3. However, only the reflection group showed lasting positive effects on attitudes (still small effect), behavioural aggression (still small effect), and compassion (still moderate effect), suggesting these effects are durable in time. Furthermore, the reflection group showed a delayed onset effect on willingness to help, whereas they were a little more willing to help in various hypothetical scenarios than the control group.
First, attitudes toward aggression (NOBAGS) only moderated one variable, blast duration. In short, while NOBAGS does not affect blast duration in the control condition, it relates to higher blast duration in the depletion condition. Although this result is theoretically consistent with the literature, it is likely a false positive given our high number of tests, the fact that this is the only variable that NOBAGS moderates, and that the p value is relatively close to 0.5. Furthermore that variable (blastduration alone) was not in the preregistration, so we might not report this finding.
Second, implicit aggression (IAT) seems to have moderated several variables. Like for NOBAGS, it also moderated blast duration, but in a three-way interaction this time. Surprisingly, for the meditation group, implicit aggression related to lower aggression, but only when depleted, whereas there was no such interaction in the waitlist group. However, as mentioned before, that variable (blastduration alone) was not in the preregistration, so we might not report this finding.
Third, implicit aggression also moderated compassionate love, again in a three-way interaction. For the waitlist group, the effect of implicit aggression clearly depends on depletion: implicit aggression relates to lower compassion in the control group (expected), but to higher compassion in the depletion group (unexpected). However, for the meditation group, the effect was absent or partly reversed.
Fourth, implicit aggression also moderated reaction time to remember an altruistic event, again in a three-way interaction. For the meditation group, higher implicit aggression relates to shorter reaction time (unexpected), unless they are depleted. However, for the waitlist group, the effect was absent or partly reversed.
In conclusion, there seems to be group differences at Time 2 and Time 3, between the experimental conditions and the control group. However, the effects in the reflection group appear not only stronger, but also more robust (i.e,. they are the only ones lasting at Time 3). Furthermore, there are also several three-way interactions between implicit attitudes, ego depletion, and group, as expected. The nature of the interactions do not seem however to perfectly align with our original predictions. A deeper exploration of the meaning of these interactions will be required.
The package references and the full script of executive code contained in this document is reproduced in the tabs below.
sessionInfo() %>% report %>% summaryThe analysis was done using the R Statistical language (v4.2.2; R Core Team, 2022) on Windows 10 x64, using the packages iterators (v1.0.14), doParallel (v1.0.17), eefAnalytics (v1.0.6), emmeans (v1.8.5), interactions (v1.1.5), sjlabelled (v1.2.0), performance (v0.10.3), see (v0.7.5), modelbased (v0.8.6), report (v0.5.7.4), foreach (v1.5.2), datawizard (v0.7.1.2), bestNormalize (v1.9.0), missForest (v1.5), rempsyc (v0.1.2.2), visdat (v0.6.0), naniar (v1.0.0), ggplot2 (v3.4.2), dplyr (v1.1.2), tidyr (v1.3.0) and psych (v2.3.3).
report::cite_packages(sessionInfo())library(rempsyc)
library(dplyr)
library(interactions)
library(performance)
library(see)
library(report)
library(datawizard)
library(modelbased)
library(ggplot2)
library(bestNormalize)
library(psych)
library(visdat)
library(missForest)
library(doParallel)
library(ggplot2)
library(emmeans)
library(sjlabelled)
library(eefAnalytics)
library(tidyr)
# Read data
data <- readRDS("Data/finaldataset_n381.rds")
report_participants(data, threshold = 1) %>% cat
# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)
sessionInfo() %>% report %>% summary
report::cite_packages(sessionInfo())
data %>%
filter(is.na(T1_Group)) %>%
nrow()
data <- data %>%
filter(!is.na(T1_Group))
report_participants(data, threshold = 1) %>% cat
# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)
data <- data %>%
mutate(part.percent = convert_na_to(part.percent, 1))
data %>%
filter(part.percent < 2/3) %>%
count(T1_Group)
data2 <- data
data <- data %>%
filter(part.percent >= 2/3)
report_participants(data, threshold = 1) %>% cat
# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)
data <- data %>%
mutate(att_check = rowSums(
select(., T1_attention1, T2_attention1, T3_attention1), na.rm = TRUE))
data %>%
count(att_check)
data <- data %>%
filter(att_check >= 2)
report_participants(data, threshold = 1) %>% cat
report_participants(data, threshold = 1, group = "T1_Group") %>% cat
# Allocation ratio
report(data$T1_Group)
report(data$T2_Condition)
get_label(data$T1_recruitment) %>% cat
report(data$T1_recruitment)
data %>%
count(T1_recruitment, sort = TRUE)
data %>%
nice_density("age", histogram = TRUE)
data %>%
count(gender, sort = TRUE)
get_label(data$T1_psycho.class) %>% cat
data %>%
count(T1_psycho.class, sort = TRUE)
get_label(data$T1_virtual.reality) %>% cat
data %>%
count(T1_virtual.reality, sort = TRUE)
get_label(data$T1_medi.experience) %>% cat
data %>%
count(T1_medi.experience, sort = TRUE, .drop = FALSE)
get_label(data$T1_disorders) %>% cat
data %>%
count(T1_disorders, sort = TRUE)
get_label(data$T1_vision) %>% cat
data %>%
count(T1_vision, sort = TRUE)
get_label(data$T1_phone) %>% cat
data %>%
count(T1_phone, sort = TRUE)
get_label(data$T1_quebec) %>% cat
data %>%
count(T1_quebec, sort = TRUE)
get_label(data$T1_student) %>% cat
data %>%
count(T1_student, sort = TRUE)
get_label(data$T1_student.program_cat) %>% cat
report(data$T1_student.program)
data %>%
count(T1_student.program_cat, sort = TRUE) %>%
filter(n > 3)
get_label(data$T1_workplace) %>% cat
report(data$T1_workplace)
data %>%
count(T1_workplace_cat, sort = TRUE) %>%
filter(n > 1)
get_label(data$T3_post.medipractice) %>% cat
data %>%
count(T3_post.medipractice, sort = TRUE)
get_label(data$T3_medipractice.which) %>% cat
report(data$T3_medipractice.which)
data %>%
count(T3_medipractice.which, sort = TRUE)
get_label(data$T3_medipractice.time) %>% cat
data %>%
count(T3_medipractice.time, sort = TRUE)
get_label(data$T3_choice.medicomp) %>% cat
data %>%
count(T3_choice.medicomp, sort = TRUE)
get_label(data$T3_followup) %>% cat
data %>%
count(T3_followup, sort = TRUE)
get_label(data$T3_consent) %>% cat
data %>%
count(T3_consent, sort = TRUE)
# Check for nice_na
nice_na(data, scales = c(
"T1_BSCS", "T1_BAQ", "T1_NOBAGS", "T1_attitude", "T1_dehumanization",
"T1_WHS", "T1_CLS", "T2_NOBAGS", "T2_attitude", "T2_dehumanization",
"T2_SMS5", "T2_PANAS", "T2_WHS", "T2_CLS", "T2_charity", "T3_NOBAGS",
"T3_attitude", "T3_dehumanization", "T3_WHS", "T3_CLS"))
# Smaller subset of data for easier inspection
data %>%
# select(manualworkerId:att_check2_raw,
# condition:condition_dum) %>%
vis_miss
# Let's use Little's MCAR test to confirm
# We have to proceed by "scale" because the function can only
# support 30 variables max at a time
library(naniar)
# We only check for the variable that had missing data, charity
# Have to divide this one in two because it is too large for the function
data %>%
select(T2_charity.moisson1_1:T2_charity.suzuki2_1) %>%
mcar_test
data %>%
select(T2_charity.conserv1_1:T2_charity.armee2_1) %>%
mcar_test
# Need logical and character variables as factors for missForest
# "Error: Can not handle categorical predictors with more than 53 categories."
new.data <- data %>%
select(-c(T1_student.program, # T1_student.program = Too many categories (> 53)
T1_already.participated, # T1_already.participated = lead to error
T3_consent)) %>% # T3_consent = lead to error
mutate(across(c(where(is.character), where(is.logical)), as.factor)) %>%
as.data.frame()
# Parallel processing
registerDoParallel(cores = 4)
# Variables
set.seed(100)
data.imp <- missForest(new.data, verbose = TRUE, parallelize = "variables")
# Total time is 2 sec (4*0.5) - 4 cores
# Extract imputed dataset
data <- data.imp$ximp
# Reverse code items 2, 4, 6, 7
data <- data %>%
mutate(across(contains("BSCS"), .names = "{col}r"))
data <- data %>%
mutate(across(ends_with(paste("BSCS", c(2, 4, 6, 7), sep = "_")),
~nice_reverse(.x, 5), .names = "{col}r"))
# Get mean BSCS
data <- data %>%
mutate(T1_BSCS = rowMeans(pick(T1_BSCS_1r:T1_BSCS_7r)))
# Get alpha & omega
data %>%
select(T1_BSCS_1r:T1_BSCS_7r) %>%
omega(nfactors = 1)
# Reverse code item 7
data <- data %>%
mutate(across(contains("BAQ"), .names = "{col}r"))
data <- data %>%
mutate(across(T1_BAQ_7, ~nice_reverse(.x, 7), .names = "{col}r"))
# Get mean BAQ
data <- data %>%
mutate(T1_BAQ = rowMeans(pick(T1_BAQ_1r:T1_BAQ_12r)))
# Get alpha & omega
data %>%
select(T1_BAQ_1r:T1_BAQ_12r) %>%
omega(nfactors = 1)
data <- data %>%
mutate(T1_attitude = rowMeans(pick(T1_attitude_1:T1_attitude_9)),
T2_attitude = rowMeans(pick(T2_attitude_1:T2_attitude_9)),
T3_attitude = rowMeans(pick(T3_attitude_1:T3_attitude_9)))
# Get alpha & omega
data %>%
select(T1_attitude_1:T1_attitude_9) %>%
omega(nfactors = 1)
data %>%
select(T2_attitude_1:T2_attitude_9) %>%
omega(nfactors = 1)
data %>%
select(T3_attitude_1:T3_attitude_9) %>%
omega(nfactors = 1)
data <- data %>%
mutate(T1_dehumanization = rowMeans(pick(T1_dehumanization_1:T1_dehumanization_9)),
T2_dehumanization = rowMeans(pick(T2_dehumanization_1:T2_dehumanization_9)),
T3_dehumanization = rowMeans(pick(T3_dehumanization_1:T3_dehumanization_9)))
# Get alpha & omega
data %>%
select(T1_dehumanization_1:T1_dehumanization_9) %>%
omega(nfactors = 1)
data %>%
select(T2_dehumanization_1:T2_dehumanization_9) %>%
omega(nfactors = 1)
data %>%
select(T3_dehumanization_1:T3_dehumanization_9) %>%
omega(nfactors = 1)
data <- data %>%
mutate(across(contains("NOBAGS"), .names = "{col}r"))
# Reverse code NOBAGS (items 1:2, 5:6, 10,12, 14:16, 20)
data <- data %>%
mutate(across(ends_with(paste0("NOBAGS.", c(
"1_1", "1_2", "3_1", "3_2", "6_1", "8_1", "10_1", "12_1", "16_1"))),
~nice_reverse(.x, 4), .names = "{col}r"))
# Get mean NOBAGS
data <- data %>%
mutate(T1_NOBAGS = rowMeans(pick(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r)),
T2_NOBAGS = rowMeans(pick(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r)),
T3_NOBAGS = rowMeans(pick(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r)))
# Get alpha & omega
data %>%
select(T1_NOBAGS.1_1r:T1_NOBAGS.16_1r) %>%
omega(nfactors = 1)
data %>%
select(T2_NOBAGS.1_1r:T2_NOBAGS.16_1r) %>%
omega(nfactors = 1)
data %>%
select(T3_NOBAGS.1_1r:T3_NOBAGS.16_1r) %>%
omega(nfactors = 1)
data <- data %>%
mutate(T1_WHS = rowMeans(pick(T1_WHS_1:T1_WHS_6)),
T2_WHS = rowMeans(pick(T2_WHS_1:T2_WHS_6)),
T3_WHS = rowMeans(pick(T3_WHS_1:T3_WHS_6)))
# Get alpha & omega
data %>%
select(T1_WHS_1:T1_WHS_6) %>%
omega(nfactors = 1)
data %>%
select(T2_WHS_1:T2_WHS_6) %>%
omega(nfactors = 1)
data %>%
select(T3_WHS_1:T3_WHS_6) %>%
omega(nfactors = 1)
data <- data %>%
mutate(T1_CLS = rowMeans(pick(T1_CLS_1:T1_CLS_21)),
T2_CLS = rowMeans(pick(T2_CLS_1:T2_CLS_21)),
T3_CLS = rowMeans(pick(T3_CLS_1:T3_CLS_21)))
# Get alpha & omega
data %>%
select(T1_CLS_1:T1_CLS_21) %>%
omega(nfactors = 1)
data %>%
select(T2_CLS_1:T2_CLS_21) %>%
omega(nfactors = 1)
data %>%
select(T3_CLS_1:T3_CLS_21) %>%
omega(nfactors = 1)
data <- data %>%
mutate(across(contains("SMS5"), .names = "{col}r"))
# Reverse code SMS5 (items 3 et 5)
data <- data %>%
mutate(across(ends_with(paste0("SMS5_", c(3, 5))),
~nice_reverse(.x, 5), .names = "{col}r"))
# Get mean SMS5
data <- data %>%
mutate(T2_SMS5 = rowMeans(pick(T2_SMS5_1r:T2_SMS5_6r)))
# Get alpha & omega
data %>%
select(T2_SMS5_1r:T2_SMS5_6r) %>%
omega(nfactors = 1)
# Get mean of PANAS
# Positive affect = 1, 3, 5, 7, 9
# Negative affect = 2, 4, 6, 8, 10
data <- data %>% mutate(
T2_PANAS_pos = rowMeans(pick(paste0("T2_PANAS_", seq(1, 9, 2)))),
T2_PANAS_neg = rowMeans(pick(paste0("T2_PANAS_", seq(2, 10, 2)))))
# Get alpha & omega
# Positive affect
data %>%
select(all_of(paste0("T2_PANAS_", seq(1, 9, 2)))) %>%
omega(nfactors = 1)
# Negative affect
data %>%
select(all_of(paste0("T2_PANAS_", seq(2, 10, 2)))) %>%
omega(nfactors = 1)
data <- data %>% mutate(
T2_Charity = rowMeans(pick(contains("charity") & ends_with("1_1"))),
T2_Familiarity = rowMeans(pick(contains("charity") & ends_with("2_1"))))
# Get alpha & omega
data %>%
select(contains("charity") & ends_with("1_1")) %>%
omega(nfactors = 1)
# Create new variable blastintensity.duration
data <- data %>%
mutate(T1_blastintensity.duration = T1_blastintensity * T1_blastduration,
T2_blastintensity.duration = T2_blastintensity * T2_blastduration,
T3_blastintensity.duration = T3_blastintensity * T3_blastduration)
data_temp <- data
data <- data2
data2 <- data
data <- data_temp
# Make list of DVs
col.list <- c("T1_blastintensity", "T1_blastduration", "T1_blastintensity.duration",
"T2_blastintensity", "T2_blastduration", "T2_blastintensity.duration",
"T3_blastintensity", "T3_blastduration", "T3_blastintensity.duration",
"T1_BSCS", "T1_BAQ", "T1_attitude", "T2_attitude", "T3_attitude",
"T1_dehumanization", "T2_dehumanization", "T3_dehumanization",
"T1_NOBAGS", "T1_WHS", "T2_WHS", "T3_WHS", "T1_CLS", "T2_CLS", "T3_CLS",
"T2_SMS5", "T2_PANAS_pos", "T2_PANAS_neg", "T2_Charity",
"T2_Familiarity", "T2_memory.altruistic", "T2_memory.aggressive")
lapply(col.list, function(x)
nice_normality(data,
variable = x,
title = x,
group = "T1_Group",
shapiro = TRUE,
histogram = TRUE))
predict_bestNormalize <- function(var) {
x <- bestNormalize(var, standardize = FALSE, allow_orderNorm = FALSE)
print(cur_column())
print(x$chosen_transform)
cat("\n")
y <- predict(x)
attr(y, "transform") <- c(attributes(y), attributes(x$chosen_transform)$class[1])
y
}
set.seed(42)
data <- data %>%
mutate(across(all_of(col.list),
predict_bestNormalize,
.names = "{.col}.t"))
col.list <- paste0(col.list, ".t")
# Group normality
named.col.list <- setNames(col.list, unlist(lapply(data, function(x) attributes(x)$transform)))
lapply(named.col.list, function(x)
nice_normality(data,
x,
"T1_Group",
shapiro = TRUE,
title = x,
histogram = TRUE))
# Plotting variance
plots(lapply(col.list, function(x) {
nice_varplot(data, x, group = "T1_Group")
}),
n_columns = 2)
plots(lapply(col.list, function(x) {
plot_outliers(data, x, group = "T1_Group", ytitle = x, binwidth = 0.3)
}),
n_columns = 2)
data %>%
as.data.frame %>%
filter(T1_Group == "Waitlist") %>%
find_mad(col.list, criteria = 3)
data %>%
as.data.frame %>%
filter(T1_Group == "Reflection") %>%
find_mad(col.list, criteria = 3)
data %>%
as.data.frame %>%
filter(T1_Group == "Meditation") %>%
find_mad(col.list, criteria = 3)
data.na <- na.omit(data[col.list])
x <- check_outliers(data.na[-c(15:16)], method = "mcd", threshold = 500)
# We have to exclude two variables that are too large following the transformations
# Otherwise we get an error:
# Error in solve.default(cov, ...) :
# system is computationally singular: reciprocal condition number = 8.99059e-20
x
plot(x)
# Winsorize variables of interest with MAD
data <- data %>%
group_by(T1_Group) %>%
mutate(across(all_of(col.list),
winsorize_mad,
.names = "{.col}.w")) %>%
ungroup()
# Update col.list
col.list <- paste0(col.list, ".w")
data <- data %>%
mutate(across(all_of(col.list), standardize, .names = "{col}.s"))
# Update col.list
col.list <- paste0(col.list, ".s")
# Let's replace original variables with the transformed variables
data[gsub(".t.w.s", "", col.list)] <- data[col.list]
# If we decide to only center variables instead of standardizing them
# (as in the preregistration), then let's do this instead
data <- data %>%
mutate(across(all_of(gsub(".t.w.s", "", col.list)),
\(x) standardize(x, center = TRUE, scale = FALSE),
.names = "{col}"))
data_temp <- data
data <- data2
col.list_temp <- col.list
col.list <- c("T1_blastintensity", "T1_blastduration", "T1_blastintensity.duration",
"T2_blastintensity", "T2_blastduration", "T2_blastintensity.duration",
"T3_blastintensity", "T3_blastduration", "T3_blastintensity.duration",
"T1_BSCS", "T1_BAQ", "T1_attitude", "T2_attitude", "T3_attitude",
"T1_dehumanization", "T2_dehumanization", "T3_dehumanization",
"T1_NOBAGS", "T1_WHS", "T2_WHS", "T3_WHS", "T1_CLS", "T2_CLS", "T3_CLS",
"T2_SMS5", "T2_PANAS_pos", "T2_PANAS_neg", "T2_Charity",
"T2_memory.altruistic", "T2_memory.aggressive")
data2 <- data
data <- data_temp
col.list <- col.list_temp
# Specify the order of factor levels for "Group".
# Otherwise R will alphabetize them.
data$T1_Group <- factor(data$T1_Group, levels = c("Meditation", "Reflection", "Waitlist"))
# Define our dependent variables
DV <- data %>% select(T2_NOBAGS:T2_Charity) %>% names
# First column (which variable)
Variable <- rep(DV, each = 3)
# Second column (which comparison)
Comparison <- rep(c("MeditationvsCTR",
"ReflectionvsCTR",
"MeditationvsReflection"),
length(DV))
# 14 == number of DV
# Make list of all formulas
formulas <- c(
"T2_NOBAGS ~ T1_Group * T1_NOBAGS",
"T2_attitude ~ T1_Group * T1_attitude",
"T2_dehumanization ~ T1_Group * T1_dehumanization",
"T2_IAT ~ T1_Group * T1_IAT",
"T2_SMS5 ~ T1_Group",
"T2_PANAS_pos ~ T1_Group",
"T2_PANAS_neg ~ T1_Group",
"T2_blastintensity ~ T1_Group * T1_blastintensity",
"T2_blastduration ~ T1_Group * T1_blastduration",
"T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
"T2_memory.altruistic ~ T1_Group",
"T2_memory.aggressive ~ T1_Group",
"T2_WHS ~ T1_Group * T1_WHS",
"T2_CLS ~ T1_Group * T1_CLS",
"T2_Charity ~ T1_Group * T2_Familiarity"
)
# Make list of all models
models.list <- sapply(formulas, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)
## Attempt with nice_lm_contrasts
x <- nice_lm_contrasts(models.list, group = "T1_Group", data = data)
x.table <- nice_table(x, highlight = TRUE)
x.table
# Make list of all formulas
formulas2 <- c(
"T3_NOBAGS ~ T1_Group * T1_NOBAGS",
"T3_attitude ~ T1_Group * T1_attitude",
"T3_dehumanization ~ T1_Group * T1_dehumanization",
"T3_IAT ~ T1_Group * T1_IAT",
"T3_blastintensity ~ T1_Group * T1_blastintensity",
"T3_blastduration ~ T1_Group * T1_blastduration",
"T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration",
"T3_WHS ~ T1_Group * T1_WHS",
"T3_CLS ~ T1_Group * T1_CLS"
)
# Make list of all models
models.list2 <- sapply(formulas2, lm, data = data, simplify = FALSE, USE.NAMES = TRUE)
## Attempt with nice_lm_contrasts
x2 <- nice_lm_contrasts(models.list2, group = "T1_Group", data = data)
x.table2 <- nice_table(x2, highlight = TRUE)
x.table2
nice_violin(data,
group = "T1_Group",
response = "T2_attitude",
obs = TRUE,
signif_annotation = c("*"),
signif_yposition = 1.5,
signif_xmin = 2,
signif_xmax = 3)
means <- estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Attitude", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(models.list$`T2_attitude ~ T1_Group * T1_attitude`)
plot(contrasts, estimate_means(models.list$`T2_attitude ~ T1_Group * T1_attitude`)) +
ylab(paste("Adjusted", "Attitude", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T2_PANAS_pos",
obs = TRUE,
comp1 = 1,
comp2 = 3,
has.d = TRUE,
d.x = 1.75,
d.y = 2)
means <- estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Positive Affect", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(models.list$`T2_PANAS_pos ~ T1_Group`)
plot(contrasts, estimate_means(models.list$`T2_PANAS_pos ~ T1_Group`)) +
ylab(paste("Adjusted", "Positive Affect", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T2_blastintensity",
obs = TRUE,
signif_annotation = c("**", "***"),
signif_yposition = c(3.5, 4),
signif_xmin = c(1, 2),
signif_xmax = c(2, 3))
means <- estimate_means(models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
models.list$`T2_blastintensity ~ T1_Group * T1_blastintensity`)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T2_blastduration",
obs = TRUE,
signif_annotation = c("**", "***"),
signif_yposition = c(3.5, 4),
signif_xmin = c(1, 2),
signif_xmax = c(2, 3))
means <- estimate_means(models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
models.list$`T2_blastduration ~ T1_Group * T1_blastduration`)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T2_blastintensity.duration",
obs = TRUE,
signif_annotation = c("*", "**"),
signif_yposition = c(3.5, 4),
signif_xmin = c(1, 2),
signif_xmax = c(2, 3))
means <- estimate_means(
models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
models.list$`T2_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T2_memory.altruistic",
obs = TRUE,
signif_annotation = c("*"),
signif_yposition = 3.5,
signif_xmin = 2,
signif_xmax = 3)
means <- estimate_means(
models.list$`T2_memory.altruistic ~ T1_Group`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list$`T2_memory.altruistic ~ T1_Group`)
plot(contrasts, estimate_means(
models.list$`T2_memory.altruistic ~ T1_Group`)) +
ylab(paste("Adjusted", "ms to remember altruistic event", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T2_CLS",
obs = TRUE,
signif_annotation = c("*", "***"),
signif_yposition = c(3.2, 2.5),
signif_xmin = c(1, 2),
signif_xmax = c(3, 3))
means <- estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Compassion", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(models.list$`T2_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list$`T2_CLS ~ T1_Group * T1_CLS`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T3_blastintensity",
obs = TRUE,
signif_annotation = c("**"),
signif_yposition = 2.5,
signif_xmin = 2,
signif_xmax = 3)
means <- estimate_means(models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)
plot(contrasts, estimate_means(
models.list2$`T3_blastintensity ~ T1_Group * T1_blastintensity`)) +
ylab(paste("Adjusted", "Blast Intensity", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T3_blastduration",
obs = TRUE,
signif_annotation = c("**"),
signif_yposition = 3,
signif_xmin = 2,
signif_xmax = 3)
means <- estimate_means(models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)
plot(contrasts, estimate_means(
models.list2$`T3_blastduration ~ T1_Group * T1_blastduration`)) +
ylab(paste("Adjusted", "Blast Duration", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T3_blastintensity.duration",
obs = TRUE,
signif_annotation = c("*"),
signif_yposition = 3,
signif_xmin = 2,
signif_xmax = 3)
means <- estimate_means(
models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)
plot(contrasts, estimate_means(
models.list2$`T3_blastintensity.duration ~ T1_Group * T1_blastintensity.duration`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T3_WHS",
obs = TRUE,
signif_annotation = c("**"),
signif_yposition = 2.5,
signif_xmin = 2,
signif_xmax = 3)
means <- estimate_means(models.list2$`T3_WHS ~ T1_Group * T1_WHS`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(
models.list2$`T3_WHS ~ T1_Group * T1_WHS`)
plot(contrasts, estimate_means(
models.list2$`T3_WHS ~ T1_Group * T1_WHS`)) +
ylab(paste("Adjusted", "Willingness to Help", "Mean")) +
theme_modern()
nice_violin(data,
group = "T1_Group",
response = "T3_CLS",
obs = TRUE,
signif_annotation = c("*", "**"),
signif_yposition = c(3.2, 2.7),
signif_xmin = c(1, 2),
signif_xmax = c(3, 3))
means <- estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)
ggplot(means, aes(x = T1_Group, y = Mean)) +
geom_line(aes(group = 1)) +
geom_pointrange(aes(color = T1_Group, ymin = CI_low, ymax = CI_high)) +
ylab(paste("Adjusted", "Compassion", "Mean")) +
theme_modern()
contrasts <- estimate_contrasts(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)
plot(contrasts, estimate_means(models.list2$`T3_CLS ~ T1_Group * T1_CLS`)) +
ylab(paste("Adjusted", "Blast Intensity * Duration", "Mean")) +
theme_modern()
data3 <- data2 %>%
mutate(part.percent = ifelse(T1_Group == "Waitlist", 0, part.percent),
part.percent = part.percent * 100,
T2_memory.altruistic = T2_memory.altruistic * -1,
T1_GroupReflection = as.numeric(T1_Group == "Reflection"),
T1_GroupMeditation = as.numeric(T1_Group == "Meditation"),
across(contains("blast"), \(x) x * -1)) %>%
as.data.frame()
caceOutput <- caceSRTBoot(
T2_attitude ~ T1_attitude + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_PANAS_pos ~ T1_GroupMeditation,
intervention = "T1_GroupMeditation",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Positive Affect (Meditation VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_blastintensity ~ T1_blastintensity + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_blastduration ~ T1_blastduration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_memory.altruistic ~ T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "ms to remember altruistic event (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_CLS ~ T1_CLS + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T2_CLS ~ T1_CLS + T1_GroupMeditation,
intervention = "T1_GroupMeditation",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_attitude ~ T1_attitude + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Attitude (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_blastintensity ~ T1_blastintensity + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_blastduration ~ T1_blastduration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Duration (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_blastintensity.duration ~ T1_blastintensity.duration + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Blast Intensity * Duration (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_WHS ~ T1_WHS + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Willingness to Help (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_CLS ~ T1_CLS + T1_GroupReflection,
intervention = "T1_GroupReflection",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Reflection VS Waitlist)")
caceOutput <- caceSRTBoot(
T3_CLS ~ T1_CLS + T1_GroupMeditation,
intervention = "T1_GroupMeditation",
compliance = "part.percent",
nBoot = 2000,
data = data3)
plot(caceOutput)
title(main = "Compassionate Love (Meditation VS Waitlist)")
# CREATE OUR DUMMY VARIABLES FOR T1_Group!
data$T1_GroupReflection <- as.numeric(data$T1_Group == "Reflection")
data$T1_GroupMeditation <- as.numeric(data$T1_Group == "Meditation")
################################################
################################################
# T2_blastintensity.duration - IAT
T2_blastintensity.duration.IAT <- lm(
T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition +
T1_GroupReflection:T2_IAT:T2_Condition + T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_blastintensity.duration.IAT)
T2_blastintensity.duration.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_blastintensity.duration.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Blast Intensity * Duration (Taylor)",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))
# T2_blastintensity - IAT
T2_blastintensity.IAT <- lm(
T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_blastintensity.IAT)
T2_blastintensity.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_blastduration - IAT
T2_blastduration.IAT <- lm(
T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
T2_IAT + T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_blastduration.IAT)
T2_blastduration.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_blastduration.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Blast Duration (Taylor)",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))
# T2_Charity - IAT
T2_Charity.IAT <- lm(
T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition + T2_Familiarity,
data = data)
check_model(T2_Charity.IAT)
T2_Charity.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Compassion - IAT
T2_CLS.IAT <- lm(
T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_CLS.IAT)
T2_CLS.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_CLS.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Compassionate Love",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))
# T2_WHS - IAT
T2_WHS.IAT <- lm(
T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT + T2_Condition +
T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_WHS.IAT)
T2_WHS.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_memory.altruistic - IAT
T2_memory.altruistic.IAT <- lm(
T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_memory.altruistic.IAT)
T2_memory.altruistic.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_memory.altruistic.IAT, pred = "T2_IAT", modx = "T2_Condition",
interval = TRUE, mod2 = "T1_GroupMeditation",
x.label = "Implicit Aggression (IAT)",
y.label = "Reaction Time (Altruistic Memory)",
mod2.labels = c("Waitlist Group", "Meditation Group"),
legend.main = "Condition") +
theme(text = element_text(size = 25))
# T2_memory.aggressive - IAT
T2_memory.aggressive.IAT <- lm(
T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_IAT +
T2_Condition + T1_GroupReflection:T2_IAT + T1_GroupMeditation:T2_IAT +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_IAT:T2_Condition + T1_GroupReflection:T2_IAT:T2_Condition +
T1_GroupMeditation:T2_IAT:T2_Condition, data = data)
check_model(T2_memory.aggressive.IAT)
T2_memory.aggressive.IAT %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_blastintensity.duration - NOBAGS
T2_blastintensity.duration.NOBAGS <- lm(
T2_blastintensity.duration ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_blastintensity.duration.NOBAGS)
T2_blastintensity.duration.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_blastintensity - NOBAGS
T2_blastintensity.NOBAGS <- lm(
T2_blastintensity ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS +
T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_blastintensity.NOBAGS)
T2_blastintensity.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_blastduration - NOBAGS
T2_blastduration.NOBAGS <- lm(
T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS +
T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_blastduration.NOBAGS)
T2_blastduration.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Make interaction plot
interact_plot(T2_blastduration.NOBAGS, pred = "T2_NOBAGS", modx = "T2_Condition",
interval = TRUE,
x.label = "Normative beliefs about aggression (NOBAGS)",
y.label = "Blast Duration (Taylor)",
legend.main = "Condition") +
theme(text = element_text(size = 25))
# T2_blastduration.NOBAGS <- lm(
# T2_blastduration ~ T1_GroupReflection + T1_GroupMeditation +
# T2_NOBAGS + T2_Condition_dum + T1_GroupReflection:T2_NOBAGS +
# T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition_dum +
# T1_GroupMeditation:T2_Condition_dum + T2_NOBAGS:T2_Condition_dum +
# T1_GroupReflection:T2_NOBAGS:T2_Condition_dum +
# T1_GroupMeditation:T2_NOBAGS:T2_Condition_dum, data = data)
# T2_blastduration.NOBAGS %>%
# nice_lm_slopes(predictor = "T2_Condition_dum",
# moderator = "T2_NOBAGS") %>%
# nice_table(highlight = TRUE)
#
# ?nice_slopes
#
# probe_interaction(T2_blastintensity.duration.NOBAGS, pred = T2_NOBAGS,
# modx = T2_Condition, digits = 3)
#
# emtrends(T2_blastintensity.duration.NOBAGS, pairwise ~ T2_Condition, var = "T2_NOBAGS")
# T2_Charity - NOBAGS
T2_Charity.NOBAGS <- lm(
T2_Charity ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition + T2_Familiarity,
data = data)
check_model(T2_Charity.NOBAGS)
T2_Charity.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# Compassion - NOBAGS
T2_CLS.NOBAGS <- lm(
T2_CLS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_CLS.NOBAGS)
T2_CLS.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_WHS - NOBAGS
T2_WHS.NOBAGS <- lm(
T2_WHS ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS + T2_Condition +
T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_WHS.NOBAGS)
T2_WHS.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_memory.altruistic - NOBAGS
T2_memory.altruistic.NOBAGS <- lm(
T2_memory.altruistic ~ T1_GroupReflection + T1_GroupMeditation +
T2_NOBAGS + T2_Condition + T1_GroupReflection:T2_NOBAGS +
T1_GroupMeditation:T2_NOBAGS + T1_GroupReflection:T2_Condition +
T1_GroupMeditation:T2_Condition + T2_NOBAGS:T2_Condition +
T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_memory.altruistic.NOBAGS)
T2_memory.altruistic.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
# T2_memory.aggressive - NOBAGS
T2_memory.aggressive.NOBAGS <- lm(
T2_memory.aggressive ~ T1_GroupReflection + T1_GroupMeditation + T2_NOBAGS +
T2_Condition + T1_GroupReflection:T2_NOBAGS + T1_GroupMeditation:T2_NOBAGS +
T1_GroupReflection:T2_Condition + T1_GroupMeditation:T2_Condition +
T2_NOBAGS:T2_Condition + T1_GroupReflection:T2_NOBAGS:T2_Condition +
T1_GroupMeditation:T2_NOBAGS:T2_Condition, data = data)
check_model(T2_memory.aggressive.NOBAGS)
T2_memory.aggressive.NOBAGS %>%
nice_lm() %>%
nice_table(highlight = TRUE)
data %>%
select(starts_with("T1_attitude")) %>%
get_label %>%
lapply(function(x) gsub(".*- ", "", x)) %>%
unlist() %>% unname
social.groups <- c("Blacks", "Homeless", "Native", "Muslims", "Refugees", "Women",
"Animals", "Elderly", "Whites")
charities <- data %>%
select(ends_with("1_1") & contains("charity")) %>%
get_label %>%
lapply(function(x) gsub(".*- ", "", x)) %>%
unlist() %>% unname
charities
regions <- c("Montreal", "Quebec", "Canada", "International")
regions
data %>%
select(T1_attitude_1:T1_attitude_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Positive Explicit Attitude",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")
data %>%
select(T1_dehumanization_1:T1_dehumanization_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Humanization",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")
data %>%
select(T2_attitude_1:T2_attitude_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Positive Explicit Attitude",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")
data %>%
select(T2_dehumanization_1:T2_dehumanization_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Humanization",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")
data %>%
select(contains("charity") & ends_with("1_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "donation") %>%
mutate(charity = factor(charity, labels = charities)) %>%
nice_violin(group = "charity",
response = "donation",
ytitle = "Amount Donated",
CIcap.width = 0.5,
obs = "jitter",
order.groups = "increasing",
xlabels.angle = 75)
data %>%
select(contains("charity") & ends_with("2_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "familiarity") %>%
mutate(charity = factor(charity, labels = charities)) %>%
nice_violin(group = "charity",
response = "familiarity",
ytitle = "Familiarity with Charity",
CIcap.width = 0.5,
obs = "jitter",
order.groups = "increasing",
xlabels.angle = 75)
data %>%
nice_scatter(predictor = "T2_Familiarity",
response = "T2_Charity",
ytitle = "Donation Amount",
xtitle = "Familiarity with Charity",
has.jitter = TRUE,
has.legend = TRUE,
has.r = TRUE,
has.p = TRUE)
data %>%
select(contains("charity") & ends_with("1_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "donation") %>%
mutate(charity = factor(charity, labels = charities),
region = case_match(
charity,
charities[1:6] ~ regions[1],
charities[7:12] ~ regions[2],
charities[13:18] ~ regions[3],
charities[19:24] ~ regions[4]
)) %>%
nice_violin(group = "region",
response = "donation",
ytitle = "Amount Donated",
obs = "jitter",
order.groups = "increasing")
data %>%
select(contains("charity") & ends_with("2_1")) %>%
pivot_longer(cols = everything(),
names_to = "charity",
values_to = "familiarity") %>%
mutate(charity = factor(charity, labels = charities),
region = case_match(
charity,
charities[1:6] ~ regions[1],
charities[7:12] ~ regions[2],
charities[13:18] ~ regions[3],
charities[19:24] ~ regions[4]
)) %>%
nice_violin(group = "region",
response = "familiarity",
ytitle = "Familiarity with Charity",
obs = "jitter",
order.groups = "increasing")
data %>%
select(T3_attitude_1:T3_attitude_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Positive Explicit Attitude",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")
data %>%
select(T3_dehumanization_1:T3_dehumanization_9) %>%
pivot_longer(cols = everything(),
names_to = "Group",
values_to = "attitude") %>%
mutate(Group = factor(Group, labels = social.groups)) %>%
nice_violin(group = "Group",
response = "attitude",
ytitle = "Humanization",
CIcap.width = 0.3,
obs = "jitter",
order.groups = "increasing")
Social Groups